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Abstract 



The main topics of this technical report are quaternions, their mathematical prop- 
erties, and how they can be used to rotate objects. We introduce quaternion math- 
ematics and discuss why quaternions are a better choice for implementing rotation 
than the well-known matrix implementations. We then treat different methods for 
interpolation between series of rotations. During this treatment we give complete 
proofs for the correctness of the important interpolation methods Slerp and Squad. 
Inspired by our treatment of the different interpolation methods we develop our own 
interpolation method called Spring based on a set of objective constraints for an 
optimal interpolation curve. This results in a set of differential equations, whose 
analytical solution meets these constraints. Unfortunately, the set of differential 
equations cannot be solved analytically. As an alternative we propose a numerical 
solution for the differential equations. The different interpolation methods are visu- 
alized and commented. Finally we provide a thorough comparison of the two most 
convincing methods (Spring and Squad). Thereby, this report provides a comprehen- 
sive treatment of quaternions, rotation with quaternions, and interpolation curves 
for series of rotations. 
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Chapter 1 



Introduction 



To animate means to "bring to life." Animation is a visual presentation of change. Traditionally 
this has been used in the entertainment business, for example Donald Duck moving in a cartoon. 
More serious applications have later been developed for physics (visualization of particle systems) 
and chemistry (displaying molecules). 

This paper treats a small part of the world of animation — animation of rotation. As a back- 
ground for the following chapters, we will in this section give an overview of how animation was 
done traditionally (i. e. before the computer), and how it is done now. The presentation is based 
on [Foley et al., 1990] and [Lasseter, 1987]. 

An animation is based on a story — a manuscript. The manuscript is used to make a storyboard, 
in which it is decided how to split the story into individual scenes. For each scene a sketch is 
made with some text describing the scene. Based on the storyboard, a series of key frames is 
produced showing the characters of the cartoon in key positions. The frames between the key 
frames can then be made from these key positions. Traditionally, the most experienced artists 
produced the key frames (and were therefore named key framers), leaving the frames in-between 
to the less experienced artists (who became known as in-betweeners). The animators produce a 
rough draft of the animation, which is presented at a pencil test. Once the draft is satisfactory, 
the final version is produced and transferred to celluloid. 

This method of animation is called key framing and has since been used in computer animation 
systems. Already in 1968 animation of 3D models was known, and the idea of using computers 
for key frame animation was used in 1971 [Burtnyk & Wein, 1971]. 

Computers are natural replacements for the in-betweeners. Given two key frames, the frames 
in-between can be generated by interpolation. Admittedly there are several problems with this 
approach: 
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• A translation between two key frames can easily be obtained by simple linear interpolation. 
When the movement consists of more key frames it is necessary to use more advanced 
curves (for example splines) to produce a smooth movement across key frames. 

• Ordinary physics cannot be used to describe how the eye perceives moving objects in a 
cartoon. Objects will change shape as they move: A ball will morph into an oval when 
bouncing fast (see figure 1.1). This will not happen automatically if a computer is used 
to animate the motion of the ball. 

• Animation of rotational movement has also been attempted using key frames and inter- 
polation. Rotation is more complex than translation, however. The problems involved in 
interpolating rotations will be treated in this paper. 




Figure 1.1: The cartoon version of a bouncing ball. 

Computer animation consists of much more than pure motion. Apart from the problems of 
interpolation of the movement there are complicated issues concerning light, sound, colors, 
camera angles, camera motion, shadows, physical properties of the objects being modelled etc. 

We limit this paper to treat methods for representation and implementation of rotation. The 
methods are mostly based on quaternions, a kind of four-dimensional complex numbers. Through 
a series of attempts to define "nice" rotation, we derive a mathematical description of rotation 
through a series of key frames. 

We will not discuss the matters mentioned in the first two bullets above or the other aspects 
mentioned (light, sound etc.) 

The main foundation for this paper is the articles [Shoemake, 1985], [Barr et al., 1992], and 
[Watt & Watt, 1992]. We do not require any knowledge of these articles. It will, however, be an 
advantage for the reader to be familiar with the common transformation methods using matrices 
and to have basic knowledge of interpolation curves in the plane (in particular splines). Some 
basic mathematics knowledge will also be advantageous (group theory, differentional calculus, 
calculus of variations and differential geometry). 
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Chapter 2 

Geometric transformations 



In this chapter we will briefly discuss selected transformations of objects in 3D. The key topic 
will be rotation but since interpolation between positions offers useful parallels to interpolation 
of rotation, we include translation. Note that these parallels serve only as inspiration for the 
rotational case — mainly because the space of translations is Euclidean while the space of 
rotations is not. This difference will be discussed in depth in the following chapters. 

2.1 Translation 

Translation is the most obvious kind of transformation: A point in space is moved from one 
position to another. Let a point P £ K 3 be denoted by a 3-tuple (x,y,z), x,y,z G M. and 
the translation by a vector (Ax, Ay, Az). Then the new position P 1 is calculated by simple 
addition: P' = (x + Ax, y + Ay, z + Az). The definition is non-ambiguous i. e. there exists only 
one translation vector that takes P to P'. 



2.2 Rotation 

Rotation in 3D is not as simple as translation and it can be defined in many ways. We have 
chosen the following definition: 

We will use the definition given by Euler's (*1707 - f 1783) theorem [Euler, 1752] — written in 
modern notation (compare with figure 2.1): 

Proposition 1. 

Let O, O' G M 3 be two orientations. Then there exists an axis I G M 3 and an angle of rotation 
9 G ] — 7r, 7r] such that O yields O' when rotated 9 about I. 

Note that the proposition states existence and does not state uniqueness. 

We will distinguish between orientations and rotations. An orientation of an object in M 3 is 
given by a normal vector. A rotation is defined by an axis and an angle of rotation. 
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Figure 2.1: Let O, O' € M 3 oe iwo orientations. Then there exists an axis I € M 3 and an angle 
of rotation 9 £ ] — tt, tt] smc/i that O yields O' when rotated 6 about I. 
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Chapter 3 



Two rotational modalities 



Euler's theorem (proposition 1) gives a simple definition of rotations. In most of the literature, 
Euler angles are used to define rotation. From these two fundamental definitions, rotation can 
be discussed mathematically in numerous ways. We will term the combination of a definition 
and a corresponding mathematical representation a rotational modality. In this report we will 
discuss the following two modalities: 

• Rotation denned by Euler angles represented by general transformation matrices. 

• Rotation denned by Euler's theorem represented by quaternions. 

The aim is of this chapter is to reach an implementation of a general rotation with each modality. 
A comparison of how the modalities implement a general rotation is given in chapter 4. Conver- 
sion between representations of rotation is discussed in appendix B. Finally, general conventions 
for rotation used in this report can be found in appendix A. 

In sections 3.1 and 3.2, Euler angles and their matrix representation are described. The descrip- 
tion is brief — the reader is assumed familiar with these topics. 

Section 3.3 gives an in-depth treatment of quaternions starting of with the basics of quaternion 
mathematics. After the introduction, it is established how quaternions can be used to represent 
rotation as denned by Euler's theorem. 

3.1 Euler angles 

The space of orientations can be parameterized by Euler angles. When Euler angles are used, 
a general orientation is written as a series of rotations about three mutually orthogonal axes in 
space. Usually the x, y, and z axes in a Cartesian coordinate system are used. The rotations 
are often called x-roll, y-roll and z-roll. 

Euler originally developed Euler angles as a tool for solving differential equations. Later Euler 
angles have become the most widely used method of parametererizing the space of orientations. 
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As we shall see below, this choice gives rise to a number of problems. If we choose to consider 
a rotation as the action performed to obtain a given orientation, Euler angles can be used to 
parameterize the space of rotations. To describe a general rotation as described in section 2.2, 
three Euler angles (#i,#2 5 #3) are required, where 8\, 62, and 63 are the rotation angles about 
the x, y, and z axes, respectively. 

The conversion from a general rotation to Euler angles is ambiguous since the same rotation 
can be obtained with different sets of Euler angles (see [Foley et al., 1990]). Furthermore, the 
resulting rotation depends on the order in which the three rolls are performed. This gives rise to 
further ambiguity but fits well with the fact that rotations in space do not generally commute 
(see appendix B). Some of the ambiguity in the conversion to Euler angles can be eliminated 
by adopting a convention of which order the rolls should be performed. In this paper, we use 
the convention described in appendix A. Introducing a convention does not, however, eliminate 
the ambiguity altogether (see chapter 4). 



3.2 Rotation matrices 



Rotation matrices are the typical choice for implementing Euler angles. For each type of roll, 
there is a corresponding rotation matrix, i. e. an x rotation matrix, a y rotation matrix, and a 
z rotation matrix. The matrices rotate by multiplying them to the position vector for a point 
in space, and the result is the position vector for the rotated point. A rotation matrix is a 
3x3 matrix, but usually homogeneous 4x4 matrices are used instead (see [Foley et al., 1990] 
for further detail). A general rotation is obtained by multiplying the three roll- matrices corre- 
sponding to the three Euler angles. The resulting matrix embodies the general rotation and can 
be applied to the points that are to be rotated. 

The three standard rotation matrices are given in homogeneous coordinates in appendix B. 

Matrix multiplication is not generally commutative. This fits well with the fact that rotations 
in space do not commute. 

Finally it should be noted that using homogeneous transformation matrices gives the only imple- 
mentation that effectively embodies all standard transformations: Translation, scaling, shearing, 
and various projection transformations. 
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3.3 Quaternions 



The second rotational modality is rotation denned by Euler's theorem and implemented with 
quaternions. Since quaternions are not nearly as well-known as transformation matrices, and 
since no good overview of the field exists, we will give a historical overview and then provide a 
thorough treatment of quaternion mathematics. 

3.3.1 Historical background 

Quaternions were invented by Sir William Rowan Hamilton (*1809 - f 1865) in 1843. Hamilton's 
aim was to generalize complex numbers to three dimensions, i. e. numbers of the form a + ifr+jc, 
where a,b, c E M and i 2 = j 2 = —1. Hamilton never succeeded in making this generalization, 
and it has later been proven that the set of three-dimensional numbers is not closed under 
multiplication. In 1966 Kenneth O. May gave the following elegant proof of this: 

Proposition 2. 

The set of three-dimensional complex numbers is not closed under multiplication. 
Proof (freely adopted from Kenneth O. May 1966): 

Assume that the usual rules of arithmetic for complex numbers hold, and that i 2 = j 2 = —1. 

The proof is by contradiction, so we assume that a closed multiplication exists. Since multi- 
plication is closed, there exist a, 6, c € R that satisfy ij = a + ib + jc. Multiplying this with 
i yields — j = — b + ia + ijc. Substituting the first equation in the second equation yields 
-j = -b + ia + (a + ib+jc)c, i. e. = (ac - b) + i(a + be) +j(c 2 + 1). Thus ac-b = 0, a + bc = 
and c 2 + 1 = 0. The equation c 2 + 1 = gives the contradiction, since c is real by assumption. 

□ 

One of Hamilton's motivations for seeking three-dimensional complex numbers was to find a 
description of rotation in space corresponding to the complex numbers, where a multiplication 
corresponds to a rotation and a scaling in the plane. 

While walking by the Royal Canal in Dublin on a Monday in October 1843, Hamilton realized 
that four numbers are needed to describe a rotation followed by a scaling. One number describes 
the size of the scaling, one the number of degrees to be rotated, and the last two numbers 
give the plane 1 in which the vector should be rotated. After this insight, Hamilton found a 
closed multiplication for four-dimensional complex numbers of the form ix + jy + k?, where 
i 2 = j 2 = k 2 = ijk = —1. Hamilton dubbed his four-dimensional complex numbers quaternions. 
The parallel to ordinary complex numbers stems from the imaginary parts. 

A quaternion is usually written [s,v],s E M, v E M 3 . Here s is called the scalar part, and 
v = (x, y, z) is the vector part. 

1 The xy plane can be rotated to any plane in xyz space through the origin by giving the rotation angles about 
the x and y axes. 
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Historical aside 

Hamilton presented quaternion mathematics at a series of lectures at the Royal Irish Academy. 
The lectures gave rise to a book [Hamilton, 1853], the full title of which (with typography as in 
the book) is: 

Lectures on Quaternions: Containing a systematic statement of 

21 mm OTatfjematical Method 

OF WHICH THE PRINCIPLES WERE COMMUNICATED IN 1843 TO THE ROYAL IRISH 
ACADEMY; AND WHICH HAS SINCE FORMED THE SUBJECT OF SUCCESSIVE COURSES 
OF LECTURES, DELIVERED IN 1848 AND SUBSEQUENT YEARS IN THE HALLS OF 
TRINITY COLLEGE, DUBLIN: WITH NUMEROUS ILLUSTRATIVE DIAGRAMS, AND WITH 
SOME GEOMETRICAL AND PHYSICAL APPLICATIONS. 

In the book (page 271) Hamilton writes (again imitating the book's typography): 

285. We know then how to interpret in two apparently different ways, which are, however, 
easily perceived to have an essential connection with each other, the following symbol of 

OPERATION, 

qQq- 1 ; 

where q may be called (as before) the operator quaternion while the symbol (suppose r) of 
the operand quaternion is conceived to occupy the place marked by the parentheses. For 
we may either consider the effect of the operation, thus symbolized, to be (as in 282, 283) 
a conical rotation of the axis of the operand round the axis of the operator, through double 
the angle thereof, in such a manner as to transport the vertex of the representative angle 
of the operand to a new position on the unit sphere, without changing the magnitude of 
that angle, nor the tensor 2 of the quaternion thus operated on: or else, at pleasure, may 
regard (by 285) the operation as causing one extremity of the representative arc of the 
same operand (r) to slide along the doubled arc of the same operator (q), without any 
change in the length of the arc so sliding, nor of its inclination to the great circle along 
which its extremity thus slides. 

The historically interested reader is referred to [Hamilton, 1899] and [Hallenberg et al., 1993] 
(only available in Danish). 

3.3.2 Basic quaternion mathematics 

In this section we will state the notation used for quaternions and establish quaternion mathe- 
matics including addition, multiplication, subtraction, and multiplication with a scalar. Finally 
we define the conjugate and the inverse of a quaternion. 

2 In modern usage = norm. 
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Notation 



We use = to mean equal by definition. Closed intervals on the real line are denoted by [a, b] = 
{ x | a < x < b,a,b,x £ M}. A semi-open interval is, for example, denoted by ]a, b] = { x \ a < 
x < b,a,b,x £ M}. The set of n times differentiable functions from A to B with continuous 
derivatives we denote C n (A, B). 

Definition 1. 

The set of quaternions is denoted H. 

Quaternions consist of a scalar part s£l and a vector part v = (x, y, z) £l 3 . We will use the 
following forms: 

Definition 2. 

Let i 2 = j 2 = k 2 = ijk = —1, ij = k and ji = — k. Then q £ H can be written: 

q = [s, v] , s £ M, v £ R 3 

= [s, (x,y,z)] , s,x,y,z GR 

= s + ix+jy + kz , s,x,y,z£M. 

We will identify the set of quaternions {[s, 0] | s £ M} with M and the set {[0, v] | v G R 3 } with 
M 3 . 

Definition 3. 

Let q,q' £ H where q = [s, (x, y, z)] and q' = [s', (x', y', z')}. The addition operator, +, is defined 
q + q' = [s, v] + [s', v'] = [s, (x, y, z)} + [s', {x' , y', z')} = (s + ix + jy + kz) + (s' + ix' +jy' + kz') 

Proposition 3. (Quaternion addition) 

Let q, q' £ H, where q = [s, v] and q' = [s', v']. Then q + q 1 = [s + s', v + v'] 

Proof of proposition 3 

q + q' = [s,v] + [s',v'] 

= (s + ix+jy + kz) + {s' + \x' +y + kz') 
= (s + s') + i(x + x')+i{y + y') + k(z + z') 
= [s + s', v + v'] 

□ 

Definition 4. 

-Lei q,q' £ H where q = s + ix +jy + kz and q' = s' + ix' + jy' + kz'. Multiplication is defined 
qq' = [s,v][s',v'] = [s, (x, y, z)][s', (x' , y' , z')] = (s + ix + jy + kz)(s' + ix' + jy' + kz') 

Proposition 4. (Multiplication of quaternions) 

Let q, q' £ H, where q = [s, v] and q' = [s', v'] . Then qq' = [ss' — v • v', v x v' + sv' + s'v], where 
■ and x denote the scalar and vector product in M 3 , respectively. 



9 



Proof of proposition 4 

From definition 2 the following identities can be obtained from simple algebra: jk = i, kj = 
— i, ik = — j and ki = j. These identities are used in: 

qq' = [s,v][s',v'] 

= (s + ix + jy + kz) (s' + ix' + jy' + kz') 
= ss' — (xx' + yy' + zz') + i(sx' + s'x + yz' — zy')+ 
j(sy' + s'y + zx 1 — xz') + k(sz' + s'z + — yx') 
= [ss' — v • v', v x v' + sv' + s'v] 

□ 

Corollary 1. (to proposition 4) 

Quaternion multiplication is not generally commutative. 

Proof of proposition 1 

We give a counter-example: ij = k, but ji = — k. 

□ 

Below we will give a number of propositions without proof. The proofs are all based on the 
principle used above: The constituent quaternions are written s + ix + jy + kz. Then, using 
simple algebra and collection of terms, the result can be written as a quaternion using definition 
2. 

We state the following properties of quaternion multiplication: 
Proposition 5. 

Let p, q,q' G H and r G M. Then: 

{pq)q' = p(qq') (Quaternion multiplication is associative.) 

p(q + q') = pq+pq' (Quaternion multiplication distributes 
(q + q')p = qp + q'p across addition.) 

Multiplying quaternions by a scalar is most easily introduced by identifying r G M with the 
quaternion [r, 0] : 

Definition 5. 

Let q G H and r G M. Multiplication by a scalar is defined 

rq = [r,0]q 

Proposition 6. (Multiplication with a scalar) 

Let q G H, where q = [s, v] and let r G M. Then rq = qr = [r, 0][s, v] = [rs, rv]. 

Note that the proposition gives that multiplication with a scalar is commutative. 

q 1 

We will use the notation - to mean —q, where q G H and r£t 

r r 

We can now introduce subtraction in the usual manner: 
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Definition 6. 

Given q, q' € H, subtraction is defined q — q' = q + ( — l)q' 
The definition gives the expected: 
Proposition 7. (Quaternion subtraction) 

Let q, q' € H, where q = [s, v] and q' = [s', v']. Then q — q' = q + (— l)q' = [s — s', v — v']. 

Corresponding to the definition of the conjugate of a complex number, we define the conjugate 
of a quaternion: 

Definition 7. 

Let q € H. Then q* is called the conjugate of q and is defined by q* = [s, v]* = [s, — v]. 
The definition gives rise to the following properties: 

Proposition 8. 

Let p,q € H. Then: 

i) (q*)* = q ii) (pq)* = q*p* Hi) (p + q)* = p* + q* iv) qq* = q*q. 

The norm of a quaternion is obtained using conjugation: 
Definition 8. 

Let p € H and let the mapping \\ ■ \\ : H r\ W be defined by \\q\\ = ^qq* . This mapping is called 
the norm and \\q\\ is the norm of q. 

That this mapping is a norm in the usual sense is shown in the corollary to proposition 9. The 
norm mapping has a number of interesting properties that are summarized in: 

Proposition 9. 

Let q,q' € H and let \\ ■ \\ : H rv M be given as is definition 8. The following equations hold: 

\\q\\ = s 2 + v • v = \J s 2 + x 2 + y 2 + z 2 (3.1) 
Ikl = \\q\\ (3-2) 
Wqq'W = hWWW (3-3) 



Proof of proposition 9 

The equations 3.1 and 3.2 can be seen directly. Equation 3.3 follows from: 

\\qq'\\ = \Jqq'{qq')* = \Jqq'q'*q* = \Jq\\q'\\ 2 q* = \fqq*\\q'\\ 2 = \f\\q\\ 2 \\q'\\ 2 = \\q\\\\q'\\ 

□ 

We will later need the inner product of two quaternions. We also want to show that the norm 
mapping is indeed a norm in the usual mathematical sense. From equation 3.1 in proposition 9 
it follows that the norm of a quaternion q can be written as it is usually obtained from the inner 
product (if q € H is identified with the corresponding vector in M 4 ). This property is formalized 
by: 
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Definition 9. 

Let q,q' E H, q = [s, v] = [s, (x, y, z)], q' = [s', v'] = [s', (x', y', z')]. The inner product is defined 
as • : H x H rv IR where q • = ss' + v ■ v' = ss' + sa;' + yy' + zz'. 

Note that the definition yields g • q = s 2 + a; 2 + y 2 + £ 2 , which gives rise to: 

Corollary 2. (to proposition 9) 

The norm of a quaternion q can be obtained by \\q\\ = ^Jq • q. Furthermore, \\ ■ \\ is a norm in 
the usual mathematical sense. 

Proof of corollary 2 

That • computes the norm squared follows directly from proposition 9 and definition 9. Now 
let q = [s, (x, y, z)] G H. If we identify q with (s, x, y, z) Gl 4 , the above method of computing 
the norm is identical to the usual Euclidean norm on IR 4 . Thus the quaternion norm is a norm 
in the usual sense. 

□ 

We will later need the following generalization of the two- and three-dimensional cases: 
Proposition 10. 

Let q,q' E H. Define q,q' as the corresponding four- dimensional vectors and let a be the angle 
between them. Then q • q' = \\q\\ \\q'\\ cos a. 

3.3.3 The algebraic properties of quaternions. 

In this section we prove that the set of quaternions H \ {[0, (0,0,0)]} is a non-Abelian group 
under quaternion multiplication. At the end of the section we give a summary of some other 
algebraic properties of quaternions. 

Definition 10. 

o 

The set of quaternions H \ {[0, (0, 0, 0)]} is written H 
We will base the discussion on the definition of a group: 
Definition 11. 

Let G be a set with an operator ■ : G x G rv G defined by (a,b) — >■ a ■ b = ab. G is a group if 

i) a(bc) = (ab)c for all a,b, c E G (The operator is associative) 

ii) Exactly one I E G exists such that la = al = (I is the neutral element) 
a for all «GG. 

Hi) For every a E G there exists an element (a~ l is the inverse element of a) 
a^ 1 E G, such that aa^ 1 = a~ 1 a = I. 

If ab = ba for all a,b E G, G is called an Abelian or commutative group. 

o 

That there exists a neutral element and inverse elements in H under quaternion multiplication 
is shown in the following two lemmas: 
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Lemma 1. 

o 

The element I = [1,0] € H is the unique neutral element under quaternion multiplication. 
Proof of lemma 1 

Let q E H be given. Proposition 6 gives ql = Iq = [Is, lv] = [s, v] = q. 

Thus 7 is a neutral element. I is also the only element that meets the requirements. To see 
this, assume that J also meets the requirements. Then IJ = 7, because J is a neutral element. 
Furthermore, I J = J, since I is a neutral element. This gives us that I = IJ = J, so I = J is 

o 

the only neutral element in H . 

□ 

Lemma 2. 

o 

Let q E H . Then there exists q^ 1 E H such that qq^ 1 = q~ 1 q = I. Furthermore q^ 1 is unique 
and given by: 




Proof of lemma 2 

o 

Let q E H be given. 

Uniqueness 

Let both pi,P2 E H be inverse to q. That p\ and p2 are equal follows from 

p 1=Pl I = p 1 (qp 2 ) = (piq)p2 = IP2 = P2 

Existence 

q* 

Let p = ,, Then 

IMr 




Thus every quaternion in H has an inverse. 

□ 

We will write | for pq~ l . Note that this is generally different from q~ l p since quaternion 
multiplication is not commutative. 

We can now state the following: 

Proposition 11. 

o 

The set His a non-Abelian group under quaternion multiplication. 
Proof of proposition 11 

Note that the set of quaternions is closed under multiplication. This follows directly from Ham- 
ilton's definition. The first requirement from the definition of a group follows from proposition 
5. The second and third requirements follow from lemmas 1 and 2. The group is not Abelian, 
since quaternion multiplication is not commutative. 

□ 
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Other algebraic properties 



The set of quaternions satisfy some other algebraic properties that are worth mentioning. These 
are given without further ado: 

• The set of quaternions is an Abelian group (H, +) under quaternion addition. 

• The set of quaternions is a non- Abelian ring (H, +,■), where + and • are quaternion 
addition and multiplication. 

3.3.4 Unit quaternions 

This section discusses a subset of the quaternion group — the set of unit quaternions. 
Definition 12. 

Let q E H . If \\q\\ = 1, then q is called a unit quaternion. We will use H\ to denote the set of 
unit quaternions. 

The set of unit quaternions constitutes a unit sphere in four-dimensional space. We shall later 
see that the set of unit quaternions play an important part in relation to general rotations. The 
following propositions lead to the important proposition 21. the following: 

Proposition 12. 

Let q = [s, v] E H\. Then there exists v' E M 3 and 6 E ] — tt, tt] such that q = [cos 6, v' sin#]. 
Proof of proposition 12 

If q = [1, 0] we let 6 = and v' can be freely chosen amongst unit vectors in M 3 . 

If q 7^ [1, 0] we let k = |v| and v' = -j:V. Then v = kv' where v' is a unit vector in IR 3 . Since q 
is a unit quaternion, we get 



The equation s 2 + k 2 = 1 describes a circle in the plane. Since a circle is also described by 
cos 2 9 + sin 2 6 = 1, there exists 6 E ] — tt, tt] such that s = cos 6 and k = sin#. All in all we get 
the desired: 



1 = 



2 



2 , 2,72/ / 2,72 

s + v ■ v = s + k v -v = s + k 



q = [ s , v] = [, 



s,v'k] = [cos 6, v' sin#] 



□ 



Two important results for unit quaternions are given in: 



Proposition 13. 

Let q,q' E H\. The following two equations hold: 

i) ll<?<?' II = 1 



ii) q 




Proof of proposition 13 

i) Wqq'W = Ikll Ik'll = 1) since \\q\ 

ii) q~ l = (/VH^I! 2 = q* , since \\q 



= \\Q 
= 1. 



1. (by equation 3.3 in proposition 9) 



□ 
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The set of unit quaternions H\ is obviously a subset of H , but definition 13 and proposition 14 

o 

give that H\ constitute a subgroup of H . 
Definition 13. 

Let G be a group and F ^ be a subset of G. F is a subgroup of G if 

i) For all a, b E F: ab E F (F is closed) 

ii) For all a £ F: a~ L £ F 



Proposition 14. 

o 

The set H\ of unit quaternions is a subgroup of the group H . 
Proof of proposition 14 

Let q, q' £ H\. Proposition 13 gives that \\qq'\\ = 1, i. e. that qq' £ Hi, and thus the first 
subgroup requirement is satisfied. Equation 3.2 in proposition 9 and proposition 13 give that 



and thereby the second subgroup requirement q 1 € H\. 

□ 



3.3.5 The exponential and logarithm functions 

We will later need quaternion versions of the real exponential and logarithm functions. The 
definitions and a few consequences of them are given here (see [Pervin & Webb, 1992] for further 
detail). 

Definition 14. 

Let q £ H\, where q = [cos#,sin#v] as in proposition 12. The logarithm function log is defined 

logg= [O,0v] 

Note that log[l, (0,0,0)] = [0, (0,0,0)] as in the real case. Note also that logg is not in general 
a unit quaternion. 

The exponential function is introduced by 
Definition 15. 

For a quaternion of the form q = [0, (9v], SeR, v£l 3 , |v| = I, the exponential function exp 
is defined by 

exp q = [cos 6, sin 6v] 

Note that the exponential and logarithm functions are mutually inverse, and that exp maps into 
Hi. 

From the above definitions we can define exponentiation for q E H\,t E M: 
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Definition 16. 

Let q E H\,t E M. Exponentiation q l is defined by 

q t = exp (Hogg) 

This gives rise to the following: 
Proposition 15. 

£ei g E -H"i, i E M. T/ien log(g*) = i log q. 
Proof of proposition 15 

log(^) = log (exp (Hogg)) = tlogq 

□ 

The following rule from M also holds for unit quaternions: 
Proposition 16. 

Let q E H\ , q = [cos 0, sin0v] and a, 6 G 1. Then 

q a q b = q a+b 

Proof of proposition 16 

q a q b = exp (a log q) exp (6 log q) 

= exp(a[O,0v])exp(6[O,0v]) 

= [cos (30, v sin a8] [cos 60, v sin 60] 

= [cos a0 cos 60 — sin a6 sin 60 ( v • v) , v cos a6 sin 60 + v cos 60 sin a6 + ( v x v) sin a0 sin 60] 

= [cos a0 cos 60 — sin a6 sin 60, v(cos a6 sin 60 + cos 60 sin a0)] 

= [cos((a + 6)0),sin((a + 6)0)v] 

= exp([0, (a + 6)0v]) 

= exp((a + 6) log(g)) 

= q a+b 

□ 

Another rule from the real numbers is (p a ) b = p ab . This rule also holds for unit quaternions: 
Proposition 17. 

Let p E Hi and a, b E M. Then (p a ) b = p ab 
Proof of proposition 17 

(p a ) b = (exp(a logp)) b = exp (6 log (exp (a logp))) = exp(6alogp) = p ab 

□ 
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One must be very careful when using exp and log as the corresponding real versions. For 
example, consider the following incorrect derivation, where p and q are unit quaternions. 

pq = exp(log(p<?)) = exp(log(p) + log(g)) = exp(log(<?) + log(p)) = exp(log(<?)) exp(log(p)) = qp 

This is inconsistent with the fact that quaternion multiplication is not commutative. The error 
lies in the second step where the rule (logpq = logp + logg) is used — this rule does not hold 
for quaternions. 

3.3.6 Rotation with quaternions 

Hamilton sought to describe rotations in space, just as complex numbers describe rotations in 
the plane. That quaternions do, in fact, perform rotation, is shown in the following propositions 
(proposition 21 in particular). 

Proposition 18. 

o 

Let p G H, p = [s, (x, y, z)] = [s, v] and let q G H . If r G M \ {0} then (rq)p(rq)^ 1 = qpq^ 1 . 
Proof of proposition 18 

Let r G M \ {0}. The inverse of rq is q~ 1 r~ 1 . Since scalar multiplication is commutative we 
can write: (rq)p(rq)~ 1 = rqpq~ 1 r~ 1 = qpq~ l rr~ l = qpq^ 1 . Thus qpq^ 1 is unchanged if q is 
multiplied by any non-zero scalar. □ 

In the propositions below, we will only consider unit quaternions, since results shown for H\ 

o 

generalize to all of H by proposition 18. 
Proposition 19. 

Let q G Hi, p = [s, v] G H. Then qpq~ l = p', where p' = [s, v'] with |v| = |v'|. 

Proof of proposition 19 

Below we write S(q) for the scalar part of q. 

The proof consists of three steps. We first show S(p') = S(p) for p G {[s, 0] | s€l} and then 
for p G {[0, v] | vgl 3 }. Finally these results are used to show the proposition for p G H. 

If p is a scalar represented as a quaternion, S(p') = S(p) follows from simple algebra. Let 
p = [s, 0], then: 

qpq- 1 = q[s, 0}q- 1 = [s, 0}qq- 1 = [s, 0] 
We have used that multiplication with a scalar commutes (proposition 6). 

Correspondingly, we will now show that the same result holds for a vector v represented as a 
quaternion [0, v]. 

The scalar part S(q) of a quaternion q can be computed by 2S(q) = q + q* . We show the 
proposition for a quaternion with in the scalar part p = [0, v]: 
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□ 



2S(qpq~ 1 ) = (qpq^ 1 ) + (qpq^ 1 )* 

= {qpq*) + (qpq*)* 

= qpq* + qp*q* (Propositions 5 and 8) 

= q(p+P*)q* (Proposition5) 

= Q(2S(p))q* 

= 2S(p) (The above result) 

= (Since p = [0, v]) 

Now let p G H, p = [s, v] = [s, 0] + [0, v]. 

qpq- 1 = q([s,0] + [0,v])q- 1 

= q[s, 0]g _1 + q[0, v]g _1 (Proposition 5 ) 
= [s, 0] + [0, v'] (The two above results) 

= [sy] 

All in all S(p') = S(p). Since q G Hi, proposition 9, equation 3.3 gives ||p'|| = 
\\q\\ \\p\\ Wq^ 1 II = Ibll- Since s is unchanged, it must be the case that |v| = |v'|. 

Corollary 3. (to proposition 19) 

Let q G Hi, p = [a, bv] G H where a, b G M and vEl 3 . // q[a, v]q* = [a, v'], i/ien <j[a, 6v]g* = 
[a,bV]. 

Proof of proposition 3 

qpq* = q[a,bv]q* 
= 9&[f,v]g* 

= 6[f,v'] (Proposition 19) 
= [aM>] 

□ 

We will later need the following useful rule: 
Proposition 20. 

Let q,p G Hi,p = [cos 0, sinflv], t G M. T/ten ijpV = (qpq*) 1 . 
Proof of proposition 20 

By corollary 3 there exists v' G M 3 such that g[cos 6, sin#v]g* = [cos 6, sin(9v']. We get 

= g(exp(i logp))g* 

= g(exp(t[0, 0\]))q* (Definition 14) 

= g(exp[O,i0v])g* 

= q([cost0,s'mtOv])q* (Definition 15) 

= [cos t6, sini#v'] (Corollary 3) 

= exp(i[O,0v']) 

= exp(i log[cos 6, sin 9v']) 

= exp(t log(qpq*)) 

= (qpq*) 1 

□ 
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We are now ready to show the main theorem of this section (inspired by [Watt & Watt, 1992]). 
Proposition 21. 

Let q E Hi, q = [cos 6, sin#n]. Let r = (x, y, z) E K 3 and p = [0, r] E H. Then p' = qpq~ l is p 
rotated 26 about the axis n. 

Proof of proposition 21 

We first show how a vector r is rotated 6 degrees about n, using sine, cosine, and the scalar 
and vector products. We then show that the same result is obtained through rotation with 
quaternions. 

Assume therefore that r is to be rotated 6 to Rr about an axis given by the unit vector n (see 
figure 3.1). 




Figure 3.1: The vector r is rotated 6 to Rr about an axis given by the unit vector n. 

The vector r can be written as a sum of two components, and rj_, where is the projection 
of r on n, and r^ is orthogonal to n (see figure 3.2). We get 

l*!! = (r • n)n, and 

rj_ = r — = r — (r ■ n)n 

To see how the rotation affects r, we place a two-dimensional coordinate system in the plane 
that is orthogonal to n and contains the points designated by r and Rr. To do this, we need a 
vector v that is orthogonal to r^ and n: 

v = n x rj_ = n x (r — (r • n)n) = nxr — n x (r ■ n)n = nxr — = nxr 
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Rr 




Figure 3.2: In the two-dimensional coordinate system orthogonal to n, (Rr)± can be written 
(Rr)± = r± cos 6 + v sin (9. 



From figure 3.2 we see that i?r's component orthogonal to n, (Rr)±, is given by 

(Rr)± = r± cos 6 + v sin 6. 

We now get: 

Rr = {Rt)\\ + {Rt) 1 _ 

= + rj_ cos 6 + v sin 9 

= (r • n)n + (r — (r • n)n) cos 6 + v sin# 

= (r • n)n — (r • n)ncos# + rcos# + vsin# 

= (1 - cos6<)(r • n)n + rcos^ + (n x r) sin6> (3.4) 

We will now examine the effect of applying a quaternion to a vector, and see that we get the 
same result as in equation 3.4. 

We now look at R q (p) = qpq^ 1 and remind the reader that p = [0, r] and that q is a unit 
quaternion [s, v]: 
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Rqip) = [sM^AbM- 1 

= [s,v][0,r][s,-v] 

= [s, v][v • r, sr — r x v] 

= [s(v • r) — v • (sr — r x v), s(sr — r x v) + (v • r)v + v x (sr — r x v)] 

= [0, s 2 r — s(r x v) + (v • r)v + v x (sr) — v x (r x v)] 

= [0, s 2 r + (v • r)v — v x (r x v) — 2s(r x v)] 

= [0, s 2 r + (v • r)v — (v • v)r + (v • r)v + 2s(v x r)] (*) 

= [0, (s 2 - v • v)r + 2(v • r)v + 2s(v x r)] 



(*) Here we use the identity vi x (v2 x V3) = (vi • V3)v2 — (vi • V2)v3 

Since q is a unit quaternion, we can write q = [cos 9, (sin#)n], where |n| = 1 (by proposition 12 
on page 14). 

Substituting this into Rqip), we get: 

R q (p) = [0, (cos 2 6'-sin 2 6'(n-n))r + 2((sin6')n-r)(sin6')n 

+2 cos #((sin#)n x r)] 
= [0, (cos 2 6»-sin 2 (9)r + (2nsin 2 <9)(n-r) 

+2cos#sin#(n x r)] 
= [0, r cos 26* + (1 - cos 26<) (n • r)n + (n x r) sin 20] 

From the above derivation, we see that the result is the same vector as in equation 3.4 except 
that the above equation has 29 instead of 9. Thus, given a unit vector n and a rotation angle 9, 
the unit quaternion [cos#,sin#n] rotates r through the angle 29 about n. □ 

As a consequence of this proposition, we get the following important corollary: 

Corollary 4. (to proposition 21) 

Any general three-dimensional rotation 9 about n, |n| = 1 can be obtained by a unit quaternion. 
Proof of corollary 4 

In the above proposition choose q such that q = [cos |, sin |n]. Thus the desired rotation is 
obtained. 

□ 

Composition of rotation is achieved by multiplying the corresponding quaternions. This is 
formalized in: 

Proposition 22. 

Let q\ , q2 £ H 1 . Rotation by q\ followed by rotation by q^ is equivalent to rotation by qiq\- 
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Proof of proposition 22 

Given p E H, the result follows directly from 



Q2{qiPQi l )q 2 1 = (<?2<?iM<?i 1 q 2 1 ) 

= (Q2Qi)p(Q*Q2) (proposition 13) 
= (<?2<?iM<?2<?i)* (proposition 8) 
= fegiM^gi) -1 (proposition 13) 

□ 

3.3.7 Geometric intuition 

We will make some observations that can aid the intuitive understanding of rotation with quater- 
nions. 

The quaternions q and q~ l 

Let q = [s, v] E H\. Then 

= q' 1 = q* = [s, -v] 

It can be useful to consider the geometric interpretation of this: The inverse of q, q~ l , rotates 
the same number of degrees as q, but the axis points in the opposite direction: 




By inverting the axis, the direction of rotation is reversed; a subsequent rotation by q cancels 
out the effect of the rotation q. 

The quaternions q and —q 

The quaternion — q represents exactly the same rotation as q (this follows from proposition 
6). This may seem surprising, but should be expected: A rotation through the angle 6 about 
the axis n can also be expressed as a rotation through the angle —6 about the axis — n. It is 
therefore aesthetically pleasing that we find both rotations on the unit quaternion sphere. The 
same duality is also found in Euler's theorem. 
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Non-unit quaternions 

It follows from proposition 18 that all quaternions on the line rq, rGt, r ^ represent the 
same rotation. 



3.3.8 Quaternions and differential calculus 

In this section we show a number of common results from differential calculus for functions 
that map into H. The results will later be used to show that certain interpolation curves are 
differentiable. 

Proposition 23. 

Let q = [cos 6, s'mOv] E Hi, t E R. Then 

j t q l = g*log(g) 

Proof of proposition 23 

The equation is shown through simple calculation of the two sides of the equation. 
The left-hand side: 

j t Q l = ^exp(tlog(g)) = ^ exp(t[0, 0v]) = ^[cosftfl), sin(*fl)v] = 0[- sm(t0), cos(^)v] 

The right-hand side: 

qHog(q) = exp(ilog(<?))log(<?) = [cos(i0), sin(i0)v][O,0v] 
= [-0sin(i0)(v • v),0cos(i0)v + 0sin(i0)(v x v)] 
= [-0sin(t0),0cos(t0)v] = 0[-sin(tfl), cos(tfl)v] 

□ 
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We also want to show the chain rule and the product rule for quaternions. We will first show 
the product rule. The purpose of this derivation is to ensure that the order of the quaternions 
in the differentiated expression is correct; it is important to make sure that this is the case since 
quaternion multiplication is not commutative. 

Proposition 24. (The product rule) 

Let f,g £ C L (M.,H). Then 



Proof of proposition 24 



-(f(t)g(t)) = lha f( x + S )9^ + S)-f(x)9(x) 

lim + + 5 )~f( x + 6 j9ix) +f{x + S)g(x) - f{x)g{x) 

<5->-0 S 

v ( ,, . n 9(x + 8) - gjx) f(x + S) - fjx) 
= hm[f(x + 6) - s + - s g(x) 

= f(x)g'(x) + f'(x)g(x) 



Proposition 25. (The chain rule) 

Let f e C l (H, H), g e C l (R, H). Then f t f{g{x)) = f'(g(x))g'(x) 
Proof of proposition 25 

We compute the derivative at an arbitrary point c£l: 

at x^c x — c 

Um (f(g(x))-f(g(c)))(g(x)-g(c))-Hg(x)-g(c)) 

x^tc X — C 

f(g(x))-f(g(c))g(x)-g(c) 



□ 



x^tc V g(x) — g(c) x — c 



f'(g(c))g'(c) 



□ 



Finally we state the following result, that has no obvious counterpart in the real numbers: 
Proposition 26. 

Let q £ C 1 (M, H\ ) , r £ C^QR, M). Since q maps into Hi, q(t) can be written [cos 6(t), v(t) sin6(t)], 
and we have 



alt 



(t) r{t) = [ -sm(r(t)e(t))(r'(t)e(t) + r(t)e'(t)y 

cos(r(t)0(t)) (r'(t)8(t) + r(t)0'(t)) v(t) + sin (r(t)0(t)) v'(t) 
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Proof of proposition 26 



— exp(r(t)log(g(t))) 

|exp(r(t)[O,v(t)0(t)]) 

|exp[O,r(i)v(i)0(i)] 
d 



dt 



[ros(r(t)0(t)),sin(r(t)0(t))v(t)] 



cos(r(t)0(t)) (r'(t)0(t) + r(t)fl'(t)) v(t) + sin (r(t)0(t)) v'(i) 



Compare this equation to the derivative of two real- valued functions u,v £ C 1 (M, M): 

—it = Ult — ti + ti log it* — u. 
tft dt v y di 

This equation is different from the equation in proposition 26, and it is unlikely that it holds 
general for quaternions. 
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3.4 An algebraic overview 



This section contains a short resume of the algebraic properties of the rotational modalities. 
The mathematical concepts are assumed to be known by the reader, and they are therefore 
not described in detail. This section elaborates on the algebra previously discussed and is only 
intended for the interested reader. 

The space of three-dimensional rotations is not a simple vector-space but a closed three-dimensional 
manifold and also a non-Abelian group (see section 3.3.3) known in the literature as 50(3) 
[Shoemake, 1994b]. Here 5 is for "special" and the 0(3) stems from the definition: 0(n) = 
{n x n matrices | 0*0 = /} — the set of orthonormal n x n matrices. 

The set of unit quaternions constitutes a subgroup of the quaternion group (see section 3.3.4). In 
the literature, H\ is also called 5 3 [Shoemake, 1985]. This subgroup constitutes a hypersphere 
in quaternion space. The spherical metric for 5 3 is equivalent to the angular metric for 50(3) 
[Shoemake, 1985]. 

Furthermore the rotation group can be projected onto the four-dimensional unit sphere of unit 
quaternions. This projection is two-to-one (see section 3.3.7): For each rotation there are 
two corresponding unit quaternions — q that is obtained directly and — q, the antipodal unit 
quaternion (see [Shoemake, 1994b], [Foley et al., 1990] and appendix B). This is because 50(3) 
has the same topology as the three-dimensional projection space called (RP 3 ), while the set 
of unit quaternions constitutes a hypersphere (5 3 ) that is topologically different from RP 3 
[Shoemake, 1994b]. 

Bearing these similarities and the small discrepancy in mind, we can see that by developing a 
smooth interpolation between unit quaternions, we get a smooth interpolation between general 
rotations. The problem is not trivial, in particular because H\ constitutes a non-Euclidean 
space, which excludes the usual interpolation methods such as splines. Our task is to find an 
equivalent interpolation curve on the surface of the four-dimensional unit sphere. 
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Chapter 4 



A comparison of quaternions, Euler 
angles and matrices 

In the previous chapter we introduced two rotational modalities: 

• Rotation denned by Euler angles represented by general transformation matrices. 

• Rotation denned by Euler's theorem represented by quaternions. 

In this chapter we will describe advantages and disadvantages of the different modalities. 

4.1 Euler angles/matrices — Disadvantages 

Traditionally homogeneous matrices have been used to represent Euler angles because the basic 
rotation matrices for rotation about the x-, y-, and z-axes are simple and well-known. This 
historically based choice has some disadvantages, though. We will discuss the disadvantages 
below. 

Lack of intuition 

Describing a general rotation as rotations about the three basis axes is not natural for an 
animator. If, for instance, the animator wants to rotate an object 30 degrees about a rotation 
axis given by the vector (1,1,1), it is quite tedious to derive the corresponding Euler angles 
about the three basis axes. 
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The order of rotation axes is important 



The user of a graphical system must express rotations in respect to a certain convention that 
defines in which order the three basis rotations are applied. Different conventions yield different 
results. For example, a rotation of (§,f,0) yields different orientations depending on which 
convention of x, y, z and y, x, z is being used. 

As an example, we examine an object in a coordinate system (see figure 4.1). By rotating the 
object in figure 4.1. i the angle f about x and then j about y the result is 4.1. ii. If the convention 
y,x,z is used instead, the rotation is done by j about y and then | about x yielding 4.1. Hi. 



Figure 4.1: In the coordinate system the x-axis points to the right, the y-axis points up and 
the z-axis point to the left. Figures ii) and Hi) show the result of applying different rotation 
conventions to the object in figure i). 

Gimbal lock 

Getting an intuitive understanding of how rotation matrices work is quite difficult. In partic- 
ular, it is difficult to predict how successive rotations about the basis axes affect each other. 
Considering that the matrix representation of Euler angles has an innate singularity in the pa- 
rameterization makes this even more difficult. It is possible to create series of rotations, where 
one degree of freedom in the rotation is lost. This situation is called gimbal lock. 

Gimbal lock is a concept originating from the air and space industry, where gyroscopes are 
used [Shoemake, 1985] [Watt & Watt, 1992] [Verplaetse, 1995]. A gyroscope basically consists 
of three concentric rings. See figure 4.2 for an illustration of this (the example is inspired by 
[McCool, 1995]). 

In 4.2.i the inner ring represents the £-axis, the center ring the y-axis and the outer ring 
represents the z-axis. A rotation about the £-axis can for example result in 4.2. ii, where the 
object is rotated approximately 45 degrees about the £-axis. If we then rotate 90 degrees about 
the y-axis, we get the situation shown in figure 4.3. 

In this situation, the x and the ^-rotation acts about the same axis. This is an example of 
gimbal lock. 




ii) 



Hi) 
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i) a) 

Figure 4.2: In the coordinate system the x-axis points up, the y-axis points to the left and the 
z-axis points to the right. The inner ring represents the x-axis, the center ring the y-axis and the 
outer ring represents the z-axis. From the starting point i) the x-axis is rotated approximately 
45 degrees in it). 




Figure 4.3: In the coordinate system the x-axis points up, the y-axis points to the left and the 
z-axis points to the right. In this situation the x and the z-rotation act about the same axis. 
This phenomenon is called gimbal lock. 

Mathematically gimbal lock corresponds to loosing a degree of freedom in the general rotation 
matrix (see appendix B): 

cos 7 sin a sin (3 — cos a sin 7 cos a cos 7 sin (3 + sin a sin 7 
cos a cos 7 + sin a sin (3 sin 7 cos a sin (3 sin 7 — cos 7 sin a 
cos (3 sin a cos a cos (3 

1 

If we let (3 = |, then a rotation with a will have the same effect as applying the same rotation 
with —7. This can be seen from the following derivation (using the addition formulas for cos 
and sin): 




R(a,(3, 7 ) 



cos (3 cos 7 
cos (3 sin 7 
— sin/3 
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-1 





cos 7 sin a — cos a sin 7 cos a cos 7 + sin a sin 7 
cos a cos 7 + sin a sin 7 cos a sin 7 — cos 7 sin a 



i?(a,§, 



7) 







-1 






1 



This expression shows that the rotation only depends on the difference a — 7 and therefore it 
has only one degree of freedom instead of two. For /3 = § changes of a and 7 result in rotations 
about the same axis. 

Implementing interpolation is difficult 

Normally the coordinates of each basis axis are interpolated independently. Thereby the inter- 
dependencies between the axes are ignored. As an example, this results in unexpected effects 
when applying simple linear interpolation. See section 6.1.1 for a treatment of this. 

Ambiguous correspondence to rotations 

Given a rotation matrix it is difficult to solve the inverse problem: What are the original 
rotations about the basis axes? In general, there is no unambiguous solution to this problem. 
See appendix A or [Shoemake & Duff, 1994] for more detail. 

In addition to this, a rotation can be represented by many different rotation matrices. All in all 
the mapping between rotations and rotation matrices is neither injective nor surjective. 

The result of composition is not apparent 

According to Euler's theorem two successive rotations can be expressed as one. The two rotation 
matrices must be composed and multiplied followed by extraction of the resulting rotation. To 
determine this rotation is tedious and in general not possible as mentioned above (see page 90 
or [Shoemake & Duff, 1994] for further detail). 

The representation is redundant 

Homogeneous matrices contain expendable information. If the matrices are to be used exclusively 
for rotation, the matrices will have zeroes for indices (4,i) and (i,4),i € {1,2,3}. In addition 
to this, the matrix uses 9 places for the 4 degrees of freedom that are necessary to describe a 
rotation according to Euler's theorem. 

On top of this numerical inaccuracies will be problematic. Since rotation matrices must be 
orthonormal there are 6 constraints (each row must be a unit vector and the columns must be 
mutually orthogonal) that must be maintained during the computations. 
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4.2 Euler angles/matrices — Advantages 



An advantage of matrix implementations is that the mathematics is well-known and that matrix 
applications are relatively easy to implement using standard packages. These advantages are 
more historically than rationally determined though. 

The main advantage of the matrix representation is the ability of the homogeneous matrix to 
represent all the other basic transformations, for example translation, scaling, projection, and 
shearing. 

4.3 Quaternions — Disadvantages 

Quaternions only represent rotation 

It is possible to implement translation using quaternions (quaternion addition can be used as 
a translation transformation interpreting the vector part of the quaternion as the translation 
vector). In [Maillot, 1990] a kind of homogeneous quaternions are denned with a multiplication 
making both translation and rotation multiplicative. 

Even though it is possible to define a homogeneous quaternion and thereby including the trans- 
lation composition, this extension is not as elegant as the homogeneous matrices. The homo- 
geneous extension appears to be ignored in the literature. Quaternions are used for rotation 
exclusively, matrices are used for all other transformations. 

Quaternion mathematics appears complicated 

Quaternions are not included in standard curriculum of modern mathematics. Some might study 
the quaternion group in algebra, but knowledge of quaternions is, in general, not widespread. 
Therefore quaternions require a bit of work in the beginning. However, quaternions should pose 
no problem for someone able to understand matrix algebra. 

4.4 Quaternions — Advantages 

Obvious geometrical interpretation 

Quaternions express rotation as a rotation angle about a rotation axis. This is a more natural 
way to perceive rotation than Euler angles. 

The obvious correspondence between Euler's theorem and rotations represented by quaternions 
gives a nice intuitive understanding of quaternions. The mapping between rotations and quater- 
nions is therefore unambiguous with the exception that every rotation can be represented by two 
quaternions. This appears to be a weakness in the quaternion representation. That q and — q 
correspond to the same rotation is on the other hand mathematically pleasing. This is because 
rotations themselves come in pairs. Given a rotation, the same rotation is obtained by rotating 
in the opposite direction about the opposite axis, (see 3.3.7 on page 22). 
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Coordinate system independency 

Quaternion rotation in not influenced by the choice of coordinate system. The user of an 
animation system does not need to worry about a certain convention of the order of rotation 
about explicit axes. 

Simple interpolation methods 

Quaternions allow elegant formulations of a range of interpolation methods. Achieving a smooth 
interpolation is therefore simpler using quaternions than Euler angles. We will give a compre- 
hensive treatment of this in chapter 6. 

Compact representation 

The representation of rotation using quaternions is compact in the sense that it is four di- 
mensional and thereby only contains the four degrees of freedom required according to Euler's 
theorem. 

In theory all non-zero quaternions can be used for rotation (by proposition 18). In practical 
applications only unit quaternions will be used. Thus, only one constraint on the representation 
must be upheld during computation compared to the six constraints on rotation matrices. 

No gimbal lock 

Since gimbal lock is innate to the matrix representation of Euler angles, this problem does not 
appear in the quaternion representation. 

Simple composition 

Rotations are easily composed when using quaternions. The composition corresponds to multi- 
plication of the involved quaternions. Rotation with q\ followed by rotation with qi is achieved 
by rotating with the quaternion qiq\- 

4.5 Conclusion 

We have stated a series of advantages and disadvantages for the two rotation modalities. Using 
Euler angles represented by matrices leads to several problems. Rotation must be expressed 
as the angles about three explicit axes, with the order being important. It is possible to en- 
counter gimbal lock and finally it is troublesome to uphold the mathematical constraints on the 
representation during calculations. 
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The quaternion representation is compact with a more natural geometrical interpretation and a 
parameterization of rotation that is not dependent on the coordinate system. 

The only real advantage of matrices is the possibility of representing all the other transforma- 
tions. 

All in all, quaternions offer the best choice for representation of rotations. 



4.6 Other modalities 

In the previous sections we have argued that rotations should be represented by quaternions 
based on Euler's theorem. However, we have only compared to one other modality — rotation 
matrices based on Euler angles. 

Obviously, other modalities are possible. For example the two modalities can be combined. It 
is possible to define rotation matrices based on Euler's theorem (an example can be found in 
[Foley et al., 1990] exercise 5.15). Thereby the problems connected to Euler angles are avoided 
yielding a better correspondence between rotation matrices and the set of rotations. There 
are still problems with this modality though. For instance, the inverse mapping from rotation 
matrices to rotations is still ambiguous. Further, the matrix representation is not well-suited for 
interpolation algorithms. For example the matrices still need to fulfill the constraints imposed 
by being orthonormal matrices. 

We will limit ourselves from discussing other modalities. This is simply because the two we 
have mentioned are by far the most important. Euler angles and matrices are the most com- 
mon modality in both the literature and in applications. We claim that quaternions are more 
appropriate. 
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Chapter 5 

Visualizing interpolation curves 



In chapter 6 we discuss a series of interpolation methods that can interpolate between two or 
more quaternions. We would like to compare these methods from a theoretical perspective. At 
the same time it is natural to compare the results of the methods in practice, to see if practice 
reflects our theoretical considerations. 

This section contains a short description of the visualization methods used and the motivation 
for each type of visualization. 

5.1 Direct visualization 

The most obvious visualization method is to apply the interpolated rotations to the object. This 
is most easily achieved by defining an object using a three-dimensional visualization tool, and 
then rotating it with the interpolated rotations. We let this visualization method produce an 
animated sequence that shows how the object is rotated. The method gives an intuitive feel 
for how the interpolation curve behaves, but it is difficult to say anything concrete about the 
smoothness of the interpolation curves or the variation in angular velocity. We therefore need 
other methods of visualization that can provide us with this information. 

5.2 Visualizing an approximation of angular velocity 

We want to visualize the angular velocity of the interpolation curve. For example it will be 
interesting to see if some of the interpolation curves have constant angular velocity. 

In the following, % denotes the «'th frame, i. e. the i'th quaternion in a discrete quaternion 
interpolation curve. 

To produce a graph of the angular velocity, we must define a function that gives an approximation 
of the angular velocity. Then all that remains is to plot this function against the interpolation 
parameter. A number of different approximations to the angular velocity can be defined based 
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in either physics or mathematics. We will base our definition in mathematics, and use that 
we have denned a norm on quaternions. We can define the distance between two quaternions 
qi,q2 to be d(q\,q2) = \\q\ — q2\\- Then the angular velocity V in the i'th quaternion % can be 
approximated by the centered average: 

v ^ = d(qi,qi-i) +d(qi,q i+ i) = \\qj - gj_ijj + \\qj - q i+ i\\ ^ ^ 

Plotting V as a function of the interpolation parameter yields a graph of an approximation to 
the angular velocity. 

We will omit the first and last key frames from the angular velocity graph, since no obvious 
angular velocity can be assigned to them. Thereby the leftmost point on the angular velocity 
graph is the angular velocity in the first interpolated frame. The remaining key frames we will 
mark with an asterisk in the velocity graph. 



5.3 Visualizing the smoothness of interpolation curves 

We would also like to visualize the interpolation curves to see, for example, how smooth they 
appear. 

Since quaternion space is four-dimensional, we cannot visualise the interpolated curves directly. 
We will always interpolate between unit quaternions, and the interpolated quaternions will 
always (with a few exceptions in chapter 6 on page 38 and 69) be unit quaternions. This 
means that we only need three dimensions to visualize the interpolation curves, because they 
lie on the surface of the unit sphere. In practice it can be difficult to visualise this space 
effectively since it must be presented via a two-dimensional media (paper or monitor). We can 
remove another dimension by interpolating between quaternions in the same three-dimensional 
hyperplane. This can be achieved by fixing one of the quaternion coordinates in all key frames. 
Thus the interpolated curves should stay inside a two-dimensional space that can be shown in 
the plane. To keep the association to the four-dimensional unit sphere, we elect to show the 
curves on the surface of the three-dimensional unit sphere (a two-dimensional sub manifold of 
M 3 ). In chapter 6 we will argue that the ideal interpolation curve will lie on the surface of the 
quaternion unit sphere. Our choice of visualization ensures that we can visually determine if 
the visualization curve stays on the surface of the unit sphere. See figure 5.3 for an example of 
this. 

The actual visualizations are produced using a ray tracer [POV, 1997]. Here we generate a large 
sphere that represents the three-dimensional unit sphere. The first key frame is shown as a 
medium-sized sphere, other key frames are shown using a bit smaller spheres, and interpolated 
points are shown using small spheres. See figures 5.1 through 5.3 below. 
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5.4 Some examples of visualization 



In this section we show some examples of the visualization of the angular velocity graphs and 
the interpolation curves. The interpolation methods are described in section 6; we only describe 
properties of the resulting visualizations here. 

The interpolation is performed on the frames given in table 5.1. The key frames are given by 
a general rotation. As noted above, we choose the rotation angle and axis such that all the 
rotations lie in the same three-dimensional hyperplane. 

In figures 5.1 through 5.3 we discuss different properties of the visualizations. Note that there 
is no obvious connection between the rotations in figure 5.1 and the points on the surface of the 
sphere. This is because the table contains general rotations, while the visualizations show the 
corresponding quaternions. Since we are only interested in the geometric shape of the curves, 
the absolute positioning of the key frames on the sphere is irrelevant. 



Rotation angle 6 E ] — tt, tt] 


Rotation axis v E M 3 


1 


(1,3,0) 


1.9 


(-1,0,0) 





(-2,1,0) 


-2 


(3,4,0) 


-1 


(-1,4,0) 


1 


(1,3,0) 



Table 5.1: Keyframes 




Figure 5.1: The interpolation curve stays on the surface of the sphere, but it is not differentiable 
in any key frames; the curve "breaks" when it passes through the key frames. The angular velocity 
graph is piecewise continuous and shows that the angular velocity is constant between keys. This 
method of interpolation is called Slerp and is described in section 6.1.5. 
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Figure 5.2: The interpolation curve is now differentiable through all key frames. Compare, for 
example, the key frame in the middle of the figure with the corresponding key frame in figure 
5.1. The angular velocity graph is continuous and assumes local minima at the key frames. This 
interpolation curve is called Squad, and it is described in section 6.2.1. 




Figure 5.3: This interpolation curve dips below the surface of the three-dimensional unit sphere. 
This means that the interpolated points are not unit quaternions, and thus the points do not lie 
on the surface of the sphere. The angular velocity graph is piecewise linear. This interpolation 
curve is called LinEuler, and is described in section 6.1.1 

In general we will illustrate the interpolation methods with the last two visualization methods: 
Sphere and graph. We have included the animated sequences for the sake of completeness, since 
it is in this context that the interpolation serves its practical purpose. 

A few examples of these animated sequences can be seen at http : //kantine . diku.dk/~myth/gif/ 
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Chapter 6 



Interpolation of rotation 



In the previous chapters we presented and discussed two rotational modalities. In this chapter 
we will investigate how well-suited the modalities are for interpolation methods. 

Gradually we will move from simple, intuitive methods to more advanced, theoretically well- 
founded interpolation methods. Parallel to the derivation of methods we will discuss what we 
perceive as criteria denning the optimal interpolation method. Our modest aim is to define and 
implement this optimal method. 



6.1 Interpolation between two rotations 

Initially we will limit ourselves to looking at interpolation between two rotations. We will use 
the rotation representations as inspiration for treating a series of simple methods. 

Each of the methods results in an interpolation curve denned as follows. Given an arbitrary 
set M we interpolate between xq E M and 27 E M parameterised by h E [0, 1]. The resulting 
interpolation curve 7 : M x M x [0, 1] rv M must satisfy the constraints: 

7(2:0, £1,0) = 2:0 
7(2:0,27,1) = 27 

6.1.1 Linear Euler interpolation: LinEuler 

The most obvious method is simply linear interpolation between two tuples of Euler angles. 
Calling this interpolation curve LinEuler, interpolation between vq = (xo,yo, zq) E K 3 and v\ = 
(27, 2/1, £1) E M 3 can be stated algorithmically using h E [0, 1] as the interpolation parameter: 

LinEuler (vo,v\,h) = vo(l — h)+vih (6.1) 
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Figure 6.1: Interpolation curve and velocity graph for Linear Euler interpolation — LinEuler. 
Between the two key frames are 300 interpolated frames. 



Figure 6.1 is not a very good illustration of the interpolation curve. The curve "lives" in the 
set of Euler angles while the illustration shows the corresponding quaternions. As described in 
chapter 5, the illustration is designed to show unit quaternions with the z coordinate equal to 
zero. The quaternions corresponding to the key frames meet these criteria but the rest of the 
interpolation curve does not. The curve consists of unit quaternions but their z coordinate is 
not generally equal to zero. Therefore the curve disappears from the surface of the sphere in 
the illustration. 

This behaviour is neither optimal 1 nor intuitively correct. 

The velocity graph shows that the animation corresponding to the interpolation will gradually 
slow down. Again, this behaviour is neither optimal nor intuitively correct. 

6.1.2 Linear Matrix interpolation: LinMat 

An alternative simple attempt is linear interpolation between rotation matrices — meaning 
linear interpolation of every single matrix element independently of the others. 

This can be stated simply algorithmically. With parameter h € [0,1], the interpolation curve 
between the rotation matrices Mq Gl 4 xK 4 and Mi Gl 4 xl 4 the curve is denned by 



As with linear Euler interpolation, the curve for linear matrix interpolation does not in general 
lie on the unit sphere, since linear interpolation between orthonormal matrices will not in gen- 
eral produce orthonormal matrices. Thus, the interpolated matrices are general homogeneous 

1 We perceive "optimal" informally so far. Later we will state a strict definition of the optimal interpolation 
curve. 



LinMat{M ,M 1 ,h) 



M (1 - ft) + Mi ft 



(6.2) 
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Figure 6.2: Interpolation curve and velocity graph for linear matrix interpolation — LinMat. 
Between the two key frames are 300 interpolated frames. 



transformation matrices containing translation, scaling, projection and other transformation el- 
ements. Thereby the interpolation can become arbitrarily wrong. For example, it is possible to 
collapse the entire object into a single point [Shoemake & Duff, 1994]. 

Our visualization methods cannot show other transformations than rotation. Therefore we are 
not able to illustrate that the interpolated matrices are not pure rotation matrices. By converting 
matrices to quaternions we preserve only the rotational part. 

In figure 6.2 we have projected the interpolation curve on to the unit quaternion sphere to show 
the pure rotational part of the interpolation curve. Finally, after these detours, we arrive at a 
quite nice interpolation curve. 

As explained, the nice illustration does not imply that the interpolation method is usable — 
only a component of the full interpolation curve is shown. The method is only discussed for 
completeness. 

6.1.3 Linear Quaternion interpolation: Lerp 

Finally, another obvious attempt is linear interpolation between rotation quaternions (called 
Lerp for linear interpolation). For qo,q\ € H and h € [0,1] this interpolation curve can be 
stated: 

Lerp(q ,qi,h) = q {l - h) + qih (6.3) 



The interpolation curve for linear interpolation between quaternions gives a straight line in 
quaternion space. The curve therefore dips below the surface of the unit sphere. Since all 
quaternions on a line through the origin give the same rotation 2 , the curve can be projected on 

2 Except the origin itself. This follows from proposition 18. 
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Figure 6.3: Interpolation curve and velocity graph for linear quaternion interpolation — Lerp. 
Between the two key frames are 300 interpolated frames. 



to the unit sphere without changing the corresponding rotations. Therefore the interpolation 
curve is normalized (figure 6.3). 

The illustration shows that Lerp is the first of the interpolation methods discussed so far to yield 
a satisfying result. Even though the interpolation curve for Lerp resembles the curve for LinMat 
(see figures 6.2 and 6.3) we must emphasize that this does not mean that the interpolations are 
alike. The illustration for LinMat is the result of a series of transformations and not a true 
image of the interpolation curve. 

Even though the interpolation curve for Lerp is nice, the velocity graph is not intuitively satis- 
fying. The speedup in the middle is due to the fact that the interpolation curve takes a "short 
cut" below the surface of the unit sphere. This is not a desired property. The intuitively correct 
velocity graph for linear interpolation is a constant function. 

6.1.4 A summary of linear interpolation 

We have attempted simple linear interpolation with the purpose of revealing which rotation 
representation is most suitable for defining interpolation curves. 

The disadvantages of LinEuler and LinMat are evident. As previously described, Euler angles 
are not the best definition for rotation and matrices are not an obvious representation. Therefore 
it is not to be expected that simple linear interpolation between pairs of Euler angles or rotation 
matrices will result in nice interpolation curves. 

In contrast the interpolation curve for Lerp is quite nice. The only problem is the varying 
velocity graph. A constant velocity is not necessarily a requirement for a curve. However, in 
this case the varying velocity is a problem since it is the result of a flaw in the method. The 
problem is that the interpolated quaternions are not unit quaternions in general. 
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For Euler angles and rotation matrices the linear interpolation can result in unacceptable in- 
terpolation curves. It does not necessarily follow from this that it is impossible to define a 
satisfying interpolation curve using these representations. It does, however, imply that it is not 
possible to define simple algorithms yielding satisfying curves. This is a direct consequence of 
the disadvantages in the Euler angle and rotation matrix representation as discussed in chapter 
4. 

On the other hand, the very simple Lerp is close to being optimal. At this stage we perceive 
the optimal interpolation curve between two key frames to be the great arc on the quaternion 
unit sphere between the two corresponding quaternions. 

Due to the previous considerations we will refrain from further attempts at deriving interpolation 
curves based on Euler angles and rotation matrices. 

6.1.5 Spherical Linear Quaternion interpolation: Slerp 

As proven in proposition 18, all quaternions on a line through the origin 3 perform the same 
rotation. However, we only want to use unit quaternions for rotation since they possess a range 
of desirable properties 4 . 

Simple linear quaternion interpolation yields a secant between the two quaternions. Therefore 
the interpolation function has larger velocity in the middle of the curve (see figure 6.3 and 
figure 6.4). Apart from this Lerp is optimal. An obvious idea is to define an interpolation 
method yielding the same interpolation curve but where the interpolated quaternions are unit 
quaternions. Instead of doing simple linear interpolation the curve should follow a great arc on 
the quaternion unit sphere from one key frame to the other. This is called great arc interpolation 
or spherical linear interpolation - Slerp. 




a) b) c) 

Figure 6.4: An illustration in the plane of the difference between Lerp and Slerp. a) The 
interpolation covers the angle v in three steps, b) Lerp — The secant across is split in four 
equal pieces. The corresponding angles are shown, c) Slerp — The angle is split in four equal 
angles. 

This interpolation can be stated algorithmically as follows ([Shoemake, 1985], [Shoemake, 1987] 
and [Shoemake, 1997]). Given go, q\ € H\ and h € [0, 1] the following four functions are equiva- 
lent expressions for spherical linear interpolation: 

3 Except the origin itself. 

4 The unit quaternion sphere is, as mentioned in chapter 3.4, equivalent to the space of general rotations. 
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Slerp(p,q,h) = p (p* q) (6.4) 

Slerp(p,q,h) = (p q*) l ~ h q (6.5) 

Slerp(p,q,h) = (q p*) h p (6.6) 

Slerp(p,q,h) = q (q* p) 1 ^ (6.7) 

Notice the pairwise symmetry yielding the intuitively correct: 

Slerp(p, q, h) = Slerp(q,p, 1 — h) 

The equivalence of the four expressions for Slerp is proven in the following proposition inspired by 
[Shoemake, 1997]. No proof of the equivalence has previously been published (this is confirmed 
by Ken Shoemake). 

Proposition 27. 

For p, q G Hi, h G M the following four expressions are equivalent: 

(1) p (p* q) h 

(2) (pq*) l - h q 

(3) (q p*) h p 

(4) q (q* p) l - h 



Proof 

First we show (1) = (3). Proposition 20 is used (page 18). 

p(p*q) h = p(p*q) h {p*p) 

= (p(p*q) h p*)p 

= (pp*qp*) h p (Proposition 20) 

= (qp*) h p 

Now we show (4) = (2): 

q(q*p) 1 ^ h = q(q*p) 1 ~ h (q*q) 
= {q{q*p) l ~ h q*)q 

= {q(q*p)q*) 1 ~ h q (Proposition 20) 
= (pq*) 1 ~ h q 

Finally we show (2) = (1) using proposition 16 from page 16 and proposition 17: 



(pq*) 1 ~ h q = (pq*){pq*)~ h q (Proposition 16) 

= (pq*){{pq*)~ 1 ) h q (Proposition 17) 

= pq*{{pq*)*) h q 

= pq*(qp*) h q 

= p(q*(qp*)q) h (Proposition 20) 

= p(q*qp*q) h 

= p{p*q) h 

□ 
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We have thus proved the equivalence of the four expressions for Slerp 5 . From now on we will 
use Slerp(p,q,h) =p(p*q) h (equation 6.4). 

That Slerp does, in fact, perform great arc interpolation on the four-dimensional quaternion 
sphere is not obvious. Often in the literature it is stated that this follows directly from the 
Lie group structure of the unit quaternions. We provide a thorough proof in proposition 28 
requiring only basic differential geometry. 

There are several different ways of proving proposition 28. One approach is to look at the 
curvature of Slerp. It is fairly easy to prove that the curvature equals one throughout the entire 
interpolation curve. Only great arcs have curvature equal one on a unit sphere. Here we use 
another approach. The key point in this proof is observing that the curve is a great arc if the 
second derivative vector is parallel (and with opposite direction) to the position vector of the 
curve, i. e. -^Slerp{p 1 q 1 h) = c Slerp(p,q,h), c < 0. This corresponds to the forces acting on 
an object describing a plane circular motion with constant angular velocity. 

Before the proof we need a lemma: 



Lemma 3. 

Let p = [s,v],qi = [s 1: (xi,yi,zi)] = [s 1: vi],q 2 
Then (jpqi) . (jpq 2 ) = \\p\\ 2 (qi • q 2 ) 

Proof of lemma 3 



[«2, {x 2 ,y 2 ,z 2 )} = [s 2 , v 2 ] € H. 



(PQi) • (PQ2) = [ss\ - v • vi, svi + siv IvxvJ • [ss 2 - V ■ v 2 , sv 2 + s 2 v + V X v 2 ] 
= S 2 SlS 2 - SSlV • v 2 - SS 2 V ■ Vi + (v ■ vi)(v • v 2 ) + 
S 2 Vl • V 2 + SS 2 V • Vl + SVl • (v X v 2 ) + 
SSlV • V 2 + SlS 2 V • V + SlV • (v X v 2 ) + 
s(v X Vl) • v 2 + s 2 (v X Vl) • v + (v X Vl) • (v X v 2 ) 
= S 2 SlS 2 + (v • vi)(v • v 2 ) + 

S 2 Vl • v 2 + SVl • (v X V 2 ) + SlS 2 V • v+ 
s(v X Vl) • v 2 + (v X Vl) • (v X v 2 ) 

: v 2 ) = (v • v)(vi • v 2 ) - (v • v 2 )(vi • v): 

S 2 S\S 2 + (v • vi)(v • v 2 ) + 
S 2 Vi • v 2 + SVl • (v X V 2 ) + SlS 2 V • v+ 
s(v x vi) • v 2 + (v • v)(vi • v 2 ) — (v • v 2 )(vi • v) 
(s 2 + v • v)sis 2 + (s 2 + v • v)vi • v 2 + 
svi • (v x v 2 ) + s(v X Vl) • v 2 

2 (<?1 • 12) + S\\ ■ (v X V 2 ) + SV 2 • (v X Vl) 



We simplify using (v x vi) 
(pqi) • (pq 2 ) 



\\P\ 



Finally, we use the identity 



v ■ (vi X v 2 ) 



x y z 

xi yi zi 
x 2 y 2 z 2 



xyiz 2 + x 2 yzi + x x y 2 z - xy 2 z x - x x yz 2 - x 2 y x z 



'Another compelling expression for Slerp is Slerp(p,q,h) 
linear interpolation p (1 — h) + q h. The equivalence with equation 6.4 can be shown as follows: p(p*q) h 



p 1 h q h . This is intuitively analogous to ordinary 



- hqh 



P 



- hqh 



This is very nice — and yet another example of how easy it is to make erroneous 



proofs with quaternions. For q,p e H and the equation (qp) h 

require commutativity. 



q p does not hold in general. This would 
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We now get: 



(PQi) • ipQ.2) = \\p\\ 2 (qi • <?2) + svi • (v x v 2 ) + sv 2 • (v x v x ) 

= IN 2 (qi - qi)+ 

s(xiyz 2 + x 2 yiz + xy 2 zi - xiy 2 z - xy\z 2 - x 2 yzi)+ 
s(x 2 yzi + xiy 2 z + xy\z 2 - x 2 y\z - xy 2 z\ - x\yz 2 ) 

= IH 2 {qi • 92) 



□ 



Proposition 28. 

The curve Slerp(p, q, h) : H\ x H\ x [0, 1] rv H\ is a great arc on the unit quaternion sphere 
between p and q. The position vector function of Slerp has constant angular velocity. 

Proof of proposition 28 

To show proposition 28 we must prove that the following four conditions are met: 

Slerp(p, q,0) = p (6.8) 

Slerp(p,q,l) = q (6.9) 

\\Slerp(p,q,h)\\ = 1, h £ [0..1] (6.10) 
d 2 

-rri7Slerp(p,q, h) = c Slerp(p,q,h), c < € M. (6-11) 
(In 1 

Conditions 6.8 and 6.9 are shown directly using the definitions for exp and log. 

Slerp(p, q,0) = p (p* q)° = p exp([0, 0]) = p[l, 0] = p 
Slerp(p,q,l) = p (p* q) 1 =p exp(log(p* q)) 
= pp* q=pp~ l q = q 

Condition 6.10 is met since exp maps into H\ (definition 15) and since the norm of a product is 
the product of the norms (proposition 9, equation 3.3): 

\\Slerp{p,q,h)\\ = \\p\\ \\{p* q) h \\ = 1 ||exp(/i log(p* q))\\ = 1 



To show condition 6.11, we need the second derivative of Slerp. Using proposition 23 we find: 

^Slerp(p,q,h) = ^p(p*q) h 

= p{p*q) h ^g{p*q) 

= Slerp(p,q,h)log(p*q) (6.12) 

d 2 

Slerp(p,q,h) = p (p* q) h \og{p* q) 2 



dh 2 



Slerp(p,q,h) log(p* q) 



2 



Condition 6.11 holds if log(p* q) 2 is a non-positive real number. Since p*, q € Hi, thenp* q € H\. 
By proposition 12 there exists flel and v G M 3 , |v| = 1 such that p* q = [cos 6, sin#v]. Then: 
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[-<9 2 v • v, <9 2 v x v] 

[-e 2 ,o] 



Thus -^Slerp{p, q,h) = c Slerp(p, q, h) where c = — 6 2 < 0. 

□ 

Having shown that Slerp(p, q,h),h € [0, 1] spans a great arc between p and q, there are still two 
possible curves depending on which direction around the unit sphere Slerp takes. The following 
proposition states that Slerp behaves as desired. 

Proposition 29. 

Letp,q € H\. Then Slerp(p,q,h), h € [0,1], spans the shortest great arc between p and q on 
the unit quaternion sphere. 

Proof of proposition 29 

Let qi = Slerp(p,q, i) and let a denote the angle between p and qi. Slerp yields the shortest 

2 2 

arc if and only if a € ] — § , § ]■ This is equivalent to cos(a) € [0, 1]. We therefore examine the 
sign of cos (a). 

Let p, q € Hi, where p = [s, v]. 

cos (a) = p • qi (Proposition 10) 

2 

= p . Slerp (p,q, 1/2) 
= p • (p (p*g)*) 



Since p*,g € -Hi it follows that p*g € -Hi- By proposition 12 there exists w G M 3 , |w| = 1 and 
ip £ ] — 7r, 7r] such that p*g = [cos(^), sin(i/>)w]. Using lemma 3 we get: 

cos a = p« [p [cos(ip), sin(-0)w] 1//2 ^ 

= p • (p exp((l/2) log[cos(-0), sin(-0)w])) 

= p . (p exp([0, (V/2)w])) 

= p • (p [cos ('0/2) , sin (-0/2) w]) 

= (p [1,0]) . (p [cos^/2),sin^/2)w]) 

= ||p|| 2 ([l,0] • [cos (-0/2) , sin (-0/2) w]) (Lemma 3) 

= ||p|| 2 cos (</>/2) 
= cos (-0/2) 



Now ip £ ] — 7r, 7r] yields cos(ip/2) > and therefore cos(a) > 0. Thus S7erp spans the shortest 
great arc between p and q. 

□ 
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We have now proven the equivalence of the four expressions for Slerp from proposition 27 
and then proven that Slerp actually produces the desired great arc. This could conclude our 
treatment of Slerp. However, the literature has traditionally avoided the use of exponentiation 
in the expression for Slerp 6 . We have encountered no problems using the expressions from 
proposition 27. However, for the sake of completeness, we will include the following expression 
for Slerp without exponentiation: 

cos (ft) = q • qi 
Q1 i M lo sm ((! ~ h)to) + gi sin(/ift) 

Slerp{q ,qi,h) = — (6.13) 

sin(ft) 

Notice that this expression is not denned for qo = ±q\. The obvious patch is Slerp(q, q, h) = q. 

The correctness of the expression above (equation 6.13) can be shown in the plane. The inter- 
polation between po and p\ (as illustrated in figure 6.5) can be written: 

/ , \ / cos(u + ht) 

q{k) = {sin(v + ht) 

The expression from equation 6.13 can — through applying the addition formulas for sin and 
cos successively — be written as: 

p sin((l - h)t) + p 1 sin(ht) 

Slerp(pQ,pi,h) = — 

sin(i) 

(cos(«) sin((l— h)t)+cos(v+t) sin(ht) \ 
iH7(t) | 
sin(«) sin((l— h)t)+sm(v+t) sin(ht) I 
sin(t) / 

(cos(«)(sin(t) cos(ht)— cos(t) sin(/it))+(cos(«) cos(t)— sin(«) sin(t)) sin(ht) \ 
sin(«)(sin(t) cos(ht)— cos(t) sin(ht)) +(s'm(v) cos(t)+cos(«) sin(t)) sin(ht) I 

/ cos(u) cos (/it) — sin(u) sin(/it) 
y sin(u) cos (ht) + cos(v) sin(/ii) 

/ cos(v + /it) 
y sin(v + ht) 

= Q{h) 

Thus, the correctness of the expression has been proven in the plane. This result can be gener- 
alized directly to four dimensions thereby proving equation 6.13. 

Slerp summarized 

The interpolation curve for Slerp (figure 6.6) forms a great arc on the quaternion unit sphere 
(as proven in proposition 28). In differential geometry terms, the great arc is a geodesic — 
corresponding to a straight line. Not only does Slerp follow a great arc (as proven in proposition 
29) it follows the shortest great arc. Thus Slerp yields the shortest possible interpolation path 
between the two quaternions on the unit sphere 7 . Furthermore Slerp has constant angular 
velocity. All in all Slerp is the optimal interpolation curve between two rotations. 



6 Since q h = exp(/ilogg) exponentiation automatically implies the use of the logarithm and exponential func- 
tions. These functions are only defined on a limited set of quaternions and they can therefore cause problems in 
conjunction with numerical inaccuracies. 

7 It should be noted that even though Slerp performs the shortest possible arc between p and q this is not nec- 



47 



a) 



b) 



Figure 6.5: Slerp in the plane, a) The interpolation goes from p to p' across the angle t. b) A 
step in the interpolation, where h £ [0, 1], q moves from p to p' . 




Angular Velocity 

Key Frames * 



150 200 250 300 
Frame nr. 



Figure 6.6: Interpolation curve and velocity graph for spherical linear quaternion interpolation 
- Slerp. Between the two key frames there are 300 interpolated frames. 



essarily optimal. Since p and —p perform the same rotation (according to equation 18 page 17), the interpolation 
between — p and q could possibly yield a shorter interpolation path. This can be established simply by comparing 
the distance between p and q, \\p — q\\, with the distance between —p and q, \\p + q\\. 
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6.2 Interpolation over a series of rotations: 
Heuristic approach 



When interpolating between two rotations Slerp is optimal. In the set of unit quaternions the 
interpolation curve of Slerp is equivalent to a straight line (the great arc). When interpolating 
between a series of rotations problems emerge: a) The curve is not smooth at the control points, 
b) The angular velocity is not constant and c) The angular velocity is not continuous at the 
controls points. 




Figure 6.7: Interpolation curve and angular velocity graph for Slerp. Between the six key frames, 
750 interpolated frames have been generated. 

A reparameterization can easily ensure continuity across the entire interpolation. Actually the 
interpolation parameter is transformed into a number of discrete frames between each pair of key 
frames. Thus a reparameterization corresponds to assigning each interval a number of frames 
relative to the size of the interval. The size of an interval can be measured as the angle 6 between 
a pair of key frames % and given by cos 6 = qi • 

Since the number of frames in each subinterval necessarily has to be an integer the angular 
velocity is, due to rounding, only approximately constant. Compare figure 6.7 and figure 6.8. 

It is not equally simple to fix the lack of smoothness. Analogously it is simple to interpolate 
between two points in the plane with a straight line, but even in the simple Euclidean space it 
is relatively complicated to create a smooth interpolation between a series of points (see figure 
6.9). 

When interpolating between a series of control points in the plane different kinds of cubic curves 
are typically used. For example this can be done with Bezier curves, which can be constructed 
quite simply. 
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Figure 6.8: Interpolation curve and angular velocity graph for Slerp. Between the six key frames, 
550 interpolated frames have been generated. The frames in the subintervals are distributed 
according to length of the interval. 



a) 



b) 



c) 



Figure 6.9: a) In the plane simple interpolation between two points is obtained by a straight 
line, b) Linear interpolation between a series of points is not differentiable in the control points. 
c) To ensure differentiability one can use cubic curves, for example splines. 




Al 



Figure 6.10: Interpolation between the points PI and P2 with a Bezier curve. The curve is 
defined as a third- order curve, where the tangent in the control points is defined by auxiliary 
points. For example the tangent in PI is defined by the auxiliary points Al and Bl (the tangent 
is Bl-Pl or Pl-Al). The differentiability is automatically assured since the curve is a third 
order curve. 
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The Bezier curve from figure 6.10 (with auxiliary points B\ and A 2 ) that interpolates between 
the control points Pi and P 2 can be expressed algorithmically (based on [Watt & Watt, 1992]) 
as three steps of linear interpolation: 



lin(xo,xi,h) = xq{1 — h) + x\h 
Bezier(P 1 ,P 2 ,B 1 ,A 2 ,h) = lin(lin(P 1 ,P 2 ,h),lin(B 1 ,A 2 ,h),2h(l-h)) 

The auxiliary points can be moved arbitrarily, which yields a change in the shape of the curve. 
The interpolation curve between PI and P2 is solely determined from the positions of auxiliary 
points Bl, A2 relative to control points PI and P2. The tangent at PI is denned by the vector 
Bl-Pl and the tangent at P2 is denned by the vector A2-P2. 

When interpolating between a series of control points, it is often desirable to ensure differentia- 
bility in the control points. This constraint can be met by making the tangents coincide in the 
control points, i.e. ensuring that Bl-Pl = Pl-Al. 

6.2.1 Spherical Spline Quaternion interpolation: Squad 

The above construction can serve as an inspiration for formulating the spherical cubic equivalent 
of a Bezier curve. This interpolation curve is called Squad (spherical and quadrangle) and was 
presented by Shoemake in [Shoemake, 1987]. 

Shoemake defines Squad as (with h € [0, 1]): 

Definition 17. 

Squad{qi,q i+ i,Si,s i+ i,h) = Slerp{Slerp{qi,q i+ i,h),Slerp{si,s i+ i,h),2h{l - h)) (6.14) 



The resulting expression for Squad is analogous to the Bezier curve, but involves spherical linear 
interpolation instead of simple linear interpolation. Bl and A2 are written Sj and Sj+i- The 
expression for Sj (equation 6.15) will be derived below. 

Correctness of Squad 

The definition of Squad is complex and therefore neither the continuity nor the differentiability 
of the resulting interpolation curve is obvious. 

Squad was originally presented in [Shoemake, 1987] which has served as general reference for a 
proof of the differentiability of Squad. [Shoemake, 1987] is no longer available 8 , and furthermore 

8 [Shoemake, 1987] is a set of course notes from SIGGRAPH 1987, and these notes are no longer available from 
University libraries or ACM. 




1 



1 



(6.15) 
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the original proof of differentiability was flawed 9 . InfKimet al., 1996] a new proof was presented. 
However, the differentiability of Squad is a consequence of a more general result in this paper 
and therefore the proof is not very thorough. In addition, the constants Sj from proposition 17 
were not derived. After corresponding with Ken Shoemake [Shoemake, 1997] we have therefore 
derived a complete proof of the differentiability of Squad. 

Proposition 30. 

Squad G C 1 

Proof 

That Squad is continuously differentiable is obvious except at the control points, since all the 
subexpressions are continuously differentiable in a given sub-interval and therefore Squad is 
continuously differentiable inside each interval. 

We must now show continuous differentiability for Squad at a given control point First 
we must show that the neighboring segments have the control points as their value at the end 
points, i. e. that Squad(qi-i,q h Si-i,s h l) = Squad{qi,q i+ i,Si,s i+ i,0): 

Squad(qi-i,qi, Sj_i, s i5 1) = Slerp{Slerp{qi-i,qi, 1), Slerp{si-i,Si, 1), 0) 

= Slerp(qi,Si,0) 
= Qi 

Squad(qi, q i+ \,s h s i+ i, 0) = Slerp(Slerp(qi, q i+1 , 0), Slerp(si, s i+ i, 0), 0) 

= Slerp(qi,Si,0) 
= Qi 

Thus Squad is continuous and has the correct value at the control points. 

We now show that Squad is continuously differentiable at a given control point. We do this by 
deriving the derivative of Squad in a given interval. Like above, we must then show that 

^Squad(qi-i,qi,Si-i,Si,l) = ^Squad(qi, q i+1 , s u s i+1 ,0) 

To find the derivative of Squad, we need the derivative of Slerp, which we get from equation 
6.12. 

We introduce the abbreviation 

gi{h) = Slerp(qi,q i+ i,h)*Slerp(si,s i+ i,h) 

Now we will find the derivative of Squad(qi, Sj, Sj+i, h) and decide how Sj and must be 
defined to ensure differentiability at the control points. 

^-Squad(qi, q i+ \,s h s i+1 ,h) = ^-Slerp(Slerp(qi, q i+ i,h), Slerp(s h s i+1 ,h),2h(l - h)) 
dt dt 

= j t (sierp{q u q l+l ,h) 9i (h) 2h ^) 

9 According to Ken Shoemake. 
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The product rule for differentiation (proposition 24) yields: 



^Squad(qi,q i+ i,Si,s i+ i,h) = ^{sierp(q h q i+ i,h)gi(h) 2h{l h ^ 

j t (Slerp( qi ,q i+1: h)^j mihf^-V + 

SIerp( qi ,q i+1 ,h) {j t i9i{h) 2h ^) 
= ^ e rp(( ?i ,( ?i+1 ,/ ^ )log(( ? *( ?i+1 ) ft (/ ^ ) 2,^ ( 1 - ,^ ) + 
Slerp(q i: q i+1: h) (J- 9i (h) 2h ^ 

Since gi{h) is a product of unit quaternions, the function values are on the unit sphere. Therefore, 
gi(h) may be written: 

gi{h) = [cos(fl 9 . (h) ),sin(fl 9 . (h) )v 9 . (h) ] 
Here ~v gi {h) 1S a un ^ vector. We can now use proposition 26 to find the derivative of gi{h) 2h ^ l ~ h ^ : 



Jt 9i{h)2h(1 ~ h) 



■ sin (2/i(l - h)0 gi(h) ) (|(2/i(l - h))0 gi(h) + 2h(l - h)±(0 giW ?) , 



+ 



cos (2/i(l - h)6 gi(h) ) (j t (2h(l - h))0 gi{h) + 2/i(l - h)-{6 gi(h) )j v 9 , 
sin (2/i(l -fr)0 9i(h) ) ^( v 9i (ft; 

■ sin (2/i(l - ((2 - 4fc)0 9i(h) + 2/i(l - h)6 gl(h) j , 

cos (2/i(l - h)e gi[h) ) ((2 - 4fc)0 9i(h) + 2h(l - h)e m ) v 9i(/i) + 
sin (2/i(l - h)0 g . {h) ^ v 9 » (h) 

Having expanded all the subexpressions of the derivative of Squad, we will now determine Sj so 
that the derivative of Squad is continuous across each control point, i. e. 

^Squad(qi-i,qi,Si-i,Si,l) = ^Squad(qi,q i+1 , s u s i+1 ,0) 

Below we write ^ {gi-i{h) 2h<yl ~ h ^) (1) for the derivative of the expression gi_i 

(hfh(l-h) applied 

to the value 1. Using algebra and rearranging, we get: 



—Squad(qi-i, q h Sj_i, s u 1) 



Slerp{qi-i,qi, 1) log (?*_!%) + 
5Zery(ft_i, gi ,l)^( gi _iW 2h ( 1 -' , ))(l) 

% iog(?r_ lft ) + ft [o, -2 (ljv^^ (i)] 

% (log (d%)- 21og([cos(0 9i _ 1 (l)),sin(0 9i _ 1 (l))v 9i _ 1 (l)])) 
ft(log(g,*_ift) -21og(^_i(l))) 
?i(log(9,*_ift) - 21og(g*Si)) 



53 



^Squad(qi, q i+1 , Sj, s i+1 , 0) = Slerp{q h q i+1 , 0) log(q*q i+1 ) + 

Slerp(q z ,q z+1 ,0)f t (g z (h) 2h ^) (0) 

= qi \og(q*q i+1 ) + ffi [0, 2 fl 9i (0)v 9i (0)] 

= % (log + 21og([cos(fl 9i (0)),sin(e 9i (0))v 9i (0)])) 

= ft(log(9?ft+i)+21og( 9i (0))) 

= 9i(log(?i ft+i) + 2 log(g*Si)) 

Thus, Sj must satisfy 

<?i(log(<?*%+i) + 2 \og(q*Si)) = %(log(g*_!%) - 2 log(g-Si)). 

Using g* = g" 1 since % € -Hi we get: 

41og(«?*Sj) = log(g*_!%) -logfe*%+i) 

'log(Ci^) " log(g?ft+i) 

exp 



% exp 



4 

log(g*_i%) - log(9,*ft+i) 



To rewrite the expression for Sj, we use the identity ((71(72)* = q^Q*- Since the constituent 

quaternions are unit quaternions, the identities q* = q -1 , and log(g*) = — log(g) also hold. 
Finally we have: 

a _ „ PYT , ( log(g?g»-i)+log(g?g»+i) 

b i — Hi ex P I ^ 

= %exp 5 l - (6.16) 



Thus Squad is continuously differentiable at the control points with si denned as above. All in 
all we have shown that Squad is continuous and continuouly differentiable across all segments. 
Further observe that the derived equation 6.16 for si is the same as equation 6.15. 

□ 



The interpolation curve generated by Squad 

The algorithmic expression for Squad yields an interpolation curve for a series of quaternions 
qo,...,qN- The expression is not defined in the first and last interval since q_\ appears in 
the expression for sq and q n +i appears in the expression for s n . Therefore it is necessary to 
define sound values for so and s n . The simplest solution is to define so = qo and sjv = qN — 
alternatively q-\ and q n +i can be defined. The choice of so and qo have little impact on the 
resulting interpolation curve and we will consider the choice an implementation detail. 

As for the interpolation curve of Slerp (figure 6.8) it is, for implementation purposes, necessary 
to produce a discrete version of the interpolation parameter and thereby selecting the number 
of interpolated frames between each key frame. For Slerp this process was simple since the arc 
length of the interpolation curve corresponds to the angle between the two involved quaternions. 
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Since the interpolation curve for Squad is rounded, it is not simple to calculate the arc length 
between each pair of key frames and thus it is not trivial to determine the number of frames 
between each pair of key frames. We choose to determine the number of frames between each 
pair of key frames relative to the distance between the the two key frames. This is not the 
optimal choice, but a simple and effective heuristic. 




Figure 6.11: Interpolation curve and angular velocity graph for Squad. Between the six key 
frames we have interpolated 550 frames. The frames have been distributed according to interval 
length. 

From figure 6.11 it is clear that Squad gives a "nice" interpolation curve. The term "nice" can 
be read as continuous and differentiable — but even this clarification is qualitatively vague: A 
continuous and differentiable curve can have any number of more or less wild twists and turns. 
From the formulation of Squad it is far from trivial to determine qualitative properties of the 
curve. Therefore we want a more objective measure from which we can define an interpolation 
curve. In this context it is no longer adequate to use qualified guesses to derive new methods 
of interpolation. However, the previously stated methods provide a good foundation for the 
development of a more general method. 

In the next sections we will seek the formulation of a more general method from a more mathe- 
matical and physical point of view. 
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6.3 Interpolation between a series of rotations: 
Mathematical approach 



So far the interpolation methods have been fairly simple and based on the rotation represen- 
tations. In principle, the interpolation is independent of which rotational modality is used to 
implement the method. Instead, the optimal interpolation curve should be denned from the 
desired properties in the space of rotations. This optimal curve can, of course, be written 
algorithmically for any sensible representation of rotation. 

The above point can be exemplified for interpolation between two rotations. The optimal inter- 
polation curve is the equivalent of a straight line in the space of rotations. This curve can be 
written algorithmically using Euler angles, but the advantage of using quaternions is that the 
curve can be stated simply — using Slerp. This is due to the previously described equivalence 
between the space of rotations and the unit quaternion sphere. 

Therefore, we will base our discussion below on the space of rotations. We will give mathemat- 
ically based demands for the optimal interpolation curve. The goal is to give an algorithmic 
description of the optimal interpolation curve. 

6.3.1 The interpolation curve 

Interpolation between rotations is denned in the space of rotations 50(3). As mentioned earlier, 
however, 50(3) and the set, Hi, of unit quaternions are topologically equivalent. We therefore 
choose to define the general interpolation in the space of unit quaternions. 

Definition 18. 

Given k control points q-i £ Hi and I = [ii,ijfc], the interpolation curve *y(t) : I rv Hi, is 
constrained by "f(ti) = q% for ti £ I. We require that ti t^, ■ ■ ■ , i ^ t^. 

6.3.2 Definitions of smoothness 

A natural requirement is that the interpolation curve is "nice." This vague term usually means 
smooth in differential geometry. However, several different definitions of smooth exist. We 
mention the following: 

Definition 19. 

Let 7(i) be the parameterization of a curve in O n (I,M n ). Smooth can then be defined in the 
following ways: 

[Madsen, 1991]: The curve >y(t) is smooth if: j(t) £ C 1 and Vi £ I : i{t) / 0. 
[Schwarz, 1989]: The curve ^(t) is smooth if: ^(t) £ C 2 . 
[Jakobsen, 1993]: The curve j(t) is smooth if: j(t) £ C°° . 

The different definitions express that it is not immediately obvious what "nice" is. We must 
therefore examine more closely which properties we want the interpolation curve to have. 
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a) b) c) d) 

Figure 6.12: Interpolation between three control points in the plane, a) Discontinuous interpo- 
lation curve, b) Continuous curve, c) C l -curve. d) C 2 -curve. 



The curve must obviously be continuous. The curve must also be differentiable. We do not want 
either "holes" or "breaks" in an animation. Thus we demand that the interpolation curve must 
be C 1 . However, it is not as obvious whether the curve should be C 2 or, for that matter, C°°. 
Figure 6.12 illustrates the different classes of curves in the plane. 

Since the control points are symmetric, it is natural to expect that the interpolation curve 
is symmetric. The illustrations in the plane clearly show that we must demand C 2 over C 1 . 
However, illustrations in the plane are not adequate ground to base this choice on, and we 
therefore postpone this decision (see section 6.3.7). 

In the first definition we find the requirement 7'(i) 7^ 0. Thus, the interpolation function is not 
allowed to contain singularities, i. e. the interpolation must not "stop." Offhand, this seems 
to be a sensible demand. However, it is possible to make sensible but contradictory demands 
to the speed of the interpolation curve. Consider, for example, animating a pendulum. The 
control points (that define the angle that the pendulum oscillates through) will lie on a straight 
line in the space of rotations. At the outer positions, it is to be expected that the pendulum 
has no velocity, corresponding to a singularity in the interpolation function. We therefore also 
postpone this decision until section 6.3.7. 

This discussion of smoothness does not bring us much further. The definitions above will not 
even allow us to differentiate between LinEuler, Lerp and Slerp. These curves are all C 2 except 
at the control points, where they are C°. Thus, it is not sufficient just to describe which class 
of functions the interpolation curve should belong to. 

6.3.3 The optimal interpolation 

Smoothness considerations do not give adequate requirements to the definition of the interpola- 
tion curve. We need an objective measure of how "nice" our curve is in rotation space. 

We will again seek inspiration in the plane (see figure 6.13). As mentioned earlier, cubic curves 
are usually used to interpolate between a series of points. There are many kinds of cubic curves; 
the Bezier curve described in section 6.2 is an example. The most common class of curves is 
splines. 
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a) b) 

Figure 6.13: Interpolating with a spline in the plane, a) Simple linear interpolation, b) Inter- 
polation with a spline. 

Discussing splines we must, as any serious project dealing with splines, write a bit about ship 
builders. Traditionally, when ship builders wanted to decide how to shape the curved parts of a 
ship, a flexible piece of metal was used. The piece of metal was fixed between a series of rivets. 
The metal piece then adapted itself to the rivets, making a nice soft curve. The ship builders 
called this tool a spline 10 . 

Viewed physically, the above description of a spline corresponds to the piece of metal achieving 
a minimum of inner tension forces subject to the constraints (the rivets). Mathematically, the 
metal piece minimizes curvature. 

In the plane, the curve can be described as follows. Given the control points (xi, yi) € K 2 , define 
7(i) € C 2 (I,M 2 ), where t is the natural parameter 11 , such that 7(f) passes through the control 
points and at the same time minimizes the expression f I \\'j"(t)\\ 2 dt. Thus the square of the 
curvature is minimized. 

This simple formulation gives a non-ambiguous definition of the interpolation curve from a 
general concept in differential geometry. We therefore choose to view the curve that minimizes 
(the square of the) curvature as the optimal interpolation curve. As we shall see below, it is not 
as simple to compute this curve in H\ as it is in the plane. 

6.3.4 Curvature in H\ 

The interpolation curve lies on H\, which is a hypersphere in quaternion space. Normally the 
curvature for a curve 7(i) is denned as ||7"(i)||, assuming that t is the natural parameter. 

The interpolation curve Slerp yields a great arc on the quaternion unit sphere. A great arc in 
Hi is equivalent to a straight line in the plane. Thus it is to expected be that a great arc does 
not have any curvature. If the curvature is computed by ||7"(i)||, the curvature will not be zero, 
but one: The curvature of the unit sphere. We therefore want to compute the curvature relative 
to the quaternion unit sphere, and not relative to quaternion space. We will call this curvature 
the local curvature. 

The definition of local curvature for a curve that lies on a surface is based in differential geometry. 
The local curvature in a point is denned as follows. Given the point on the surface, a coordinate 
system (a map) is placed in the tangent plane. The local curvature of the curve is now the 
curvature of the curve projected onto the tangent plane. This is also called tangential curvature. 

10 The ambitious project will obviously also note that modern-day architects use a refined version of the ships- 
builder's spline to draw curves. According to an architect, however, this is a myth. 
n i. e. the parameter defined from the curve length. 
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In differential geometry a great arc (the "straight line" ) is called a geodesic. Projected onto the 
tangent plane, the geodesic becomes a straight line. Thus we see, as expected, that a great arc 
does not have any local curvature. 




a) b) 

Figure 6.14: The division of^"{t) — an analogy in the plane, a) The position vector for the 
curve 7(i) lies on the surface of the quaternion unit sphere. The tangential plane is orthogonal 
to the position vector, b) 7"(i) can be split in a component 7"(i) (in the tangential plane) and 
a component 7" (t) (parallel with the position vector). 



If a curve 7(f) lies on the surface of Hi, we can split 7"(i) into two parts (see figure 6.14): a 
component, 7"(i), in the tangential plane , and a component, 7" (t), orthogonal to the tangential 
plane. The desired part of the curvature is 7"(i). Thus, the local curvature k of the curve 7 can 
be obtained: 

K(%t) = H>(t)\\ = \\Y(t)-^(t)\\ 

Since 7(f) lies on the surface of the unit sphere, 7" (t) will be parallel to 7(i). Thus we can find 
7o(i) by projecting 7"(i) onto 7(i). Thus: 

7(j l 7(j hwii; hwii 

= l"(t) - ( 7 "(i) . 7 (t))7(t) 

Definition 20. 

Given j(t) € C 2 (I, i/ie /oca/ curvature k(j, t) is defined: 

Kh,t) = \\Y(t)-(Y(t)-i(t)h(t)\\ 

Note that t is not necessarily the natural parameter. Thus the above definition is not correct 
in a strict differential geometric sense 12 . However, we are not only interested in the shape 
(curvature) of the interpolation curve. We also want a "nice" angular velocity function, i. e. one 
that minimizes angular acceleration (corresponding, from a physical viewpoint, to minimizing 
the energy). In the above expression the angular acceleration is automatically included, exactly 
because we do not reparameterize to the natural parameter. 



12 The curvature for a curve parameterised on the natural parameter can be written k(j, t) = H7" (t)||. In general 
for a smooth curve, the correct expression from differential geometry is: 

7 "(t) i'(t)(i"(t) • Y(t)) 



k(7,*) 



IIYWII 2 b'WII 4 
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6.3.5 Minimizing curvature in Hi'. Continuous, analytical solution 



We defined the optimal interpolation curve in Hi as the curve that minimizes the square of 
the curvature, but with the restriction that it must pass through the control points. We then 
defined the relevant expression for curvature in H\. We thus get the following formulation of 
the problem: 

Given the control points qi,...,q^ € Hi we seek 7 (t) € C k (I,Hi) such that there exist 
ti, . . . ,ttf £ (I) that satisfy 7 (i$) = and such that this expression is minimized: 

K{ 1 )= / ||k( 7 ,/ 1 )|| 2 ^ (6.17) 

The problem of minimizing an integral of a function is called a calculus of variations problem. 
Below, we will outline the basic method used for solving problems in the calculus of variations. 

A necessary condition for if (7) to attain a minimum is that if' (7) = 0. For J £ I and 
tp £ C k (I,Hi) we look at the "derivative": 

^*h±m^m =0 ( ,. lg) 

5-s-O S 

In the above expression, Sip is called the variation of 7, and the function 7 + Sip is called the 
comparison function. A function 7 (/i) £ C k (I,Hi) that satisfies the boundary conditions (i. e. 
liU) = Qi for i = 1, . . . , N) is called an admissible function. 

We will now derive the requirements a solution to the variation problem must meet. We therefore 
assume that both the solution 7 and the comparison function 7 + Sip are admissible, and that 
if (7) is minimal. 

For the comparison function 7 + Sip to be admissible, it must be the case that 

V>(*i) = V>(*2) = - =rl>{t N ) = (6.19) 

Furthermore, since 7 and 7 + <5i/> lie on the surface of the unit sphere, we have that ||7|| 2 = 1 
and H7 + 8ip\\ 2 = 1. Thus we have 13 : 



|2 



ii 7 r + <Jiivir+2<j(7-V') 

l + ^(^||V|| 2 + 27 . -0) 



7 • -0 = -^H^ll 2 ( 6 - 20 ) 



13 In the derivations below, it is possible to ignore the fact that the constituent expressions contain quaternion 
functions. This is due to the fact that quaternion multiplication is not used (only the scalar product, • , is used). 
Therefore, there are no problems with commutativity. 
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We want to use our knowledge of the "derivative" and examine if (7 + Sip) — if (7): 

if (7 + Sip) - if (7) 

= / N \\k( 7 + SiP, h)\\ 2 -||«( 7 , Mil 2 dh 
Jti 

= f N 11(7 + SiP)" - ((7 + Sip)" . (7 + SiP))(l + m\ 2 - - (V • 7)7|| 2 dh 

= / IIV + Sip" - (7" • 7 + <Jt" • ip + <*7 • V" + ^ V • -0) (7 + Sip) || 2 
- Il7" " (7" • 7)7|| 2 dh 



Jti 



7" - (7" • 7)7] + W - (7" • 7)V> - (7" • </>h " (</>" • 7)7] + S 2 [...]+ S% . .] 

It" - (V • 7)tII 2 ^ 



In the above expression, we have collected terms according to the exponent of S. The goal of 
the above derivations is to find the "derivative." In the expression for if (7 + Sip) — if (7) terms 
multiplied by S 2 and 5 3 can be removed because they disappear when examining the limit, after 
being divided by S. Thus the expression can be written: 



ftN 

K(-f + Sip)-K(-f) = / \\A + SB\\ 2 - \\A\\ 2 dh 

Jh 



h 

/ f2||„n2 



S l \\B\\ z + 2SA . B dh 

h 



Here A = 7" — (7" • 7)7 and B = ip" — (7" • *y)ip — (7" • ip)^ — (ip" • 7)7. Again we may omit 
terms multiplied by S 2 , and we can rewrite the expression: 



if (7 + Sip) - if (7) = / 2<f [7" - (7" • 7)7] • [1P" - (7" • j)ip - (7" • iPh - (7 • ip"h] 

ftN 

7 7" • r - (7" • 7) (V • </>) - (7" • vo(t" • 7) - (7 • ^")(t" • 7) 



, '/// 

2S f ' 
'ti 

- (7" • 7) (7 • </>") + (7" • 7) (7" • 7) (7 • rl>) + (7" • 7) (7" • V0(7 ' 7) 
+ (7" • 7) (7 • ^")(7 • 7) d/i 



We can now use the fact that 7.7= || 7 || 2 = 1, since 7 lies on the surface of the unit sphere. 

■II 

"ijv 



Furthermore, equation 6.20 is used to rewrite 7 • ip = — ^\\ip\\ 2 - We then get: 



ftN X 

if (7 + SiP) - if (7) = 25/ 7" • V>"-(7" -7)(7" • V> + 7- </>") + (V ' 7) 2 (-dl^H 2 ) ^ 

= 25 / 7" . -0" - (7" • 7) (7" • ip + 7 • -0") dh 
Jh 
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We have yet again ignored terms multiplied by 5 2 . 

After numerous rewritings, we have isolated 5 as a single factor. This factor will disappear in 
the expression for the derivative when divided by 5. We therefore now want to isolate terms 
containing This is attempted using partial integration. The expression is rewritten as follows: 

N-l 

K( 7 + 5^)-K( 7 ) = 25 L i 

i=l 

Li = 1" • V>"-(7" •7)(V / • </> + 7 •</>") dh 

= (7" " (7" • 7)7) • r ~ (7" • 7) (7" • </>) dh 

= [(7" - ii" • 7)7) • 

- /^(J^" " (7" • 7)7}) • </>' - (7" • 7) (7" • </>) dh 

J ti 

= [(7" - (V • 7)7) • - k J-w - (7" • 7)7}) • 



+ J ( J^{7" " (7" • 7)7}) • </> " (7" • 7) (7" ' </>) ^ 



Note that the above expression requires that 7 is four time s differentiable. Now we can use that 
ifi is zero in all the control points (equation 6.19) to see that the second term is zero: 

Li = [(7" - (7" • 7)7) • </>t +1 + J^{7" - (7" • 7)7}) • </> - (7" • 7) (7" • </>) dfc 

We again consider the whole expression: 

N-l 

K{j + St/)) - K{j) = 25 ^ Li 

i=i 

N-l 

= 25 [(7" - (V • 7)7) • 

2 = 1 

' r !< t-!l \ ./. /.,// ..\l-.ll 



+2S E / (^2 {V' - (7" • 7)7}) • </> - (/ • 7) (7" • </>) dfc 

= 2«j[( 7 " - ( 7 " • 7)7) • 

N-l rU+1 d 2 

+2S E / (^2 " (7" ' 7)7} " (7" • 7)7") • </> dfc 

The last rewriting uses continuity in the control points. We can find the "derivative" 
K (l + m K (l ) = 2lil „_ h „. lh) .^ 

d^O 

+ 2 E / (^2 ^" - (V • 7)7} - (7" • 7)7") • </> dfc 

2 = 1 
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Since the derivative must be zero (equation 6.18) for any ifi, we get the following requirements 
to the solution: 



7 € <? 4 ((/),#i) 

= 7"(ti) - h"(h) . 7 (fi))7(*i) 

= 7"(tjv) - (V'(iiv) • l(t N )h(t N ) 

(l" -7)7" = J^{t" " (7" • 7)7} 

= 7"" " (7" • 7)"7 " 2(7" • 7)Y " (V • 7)7" 

The second and third requirements are equivalent with the local curvature being zero at the 
end-points. The last requirement can be rewritten using 7.7= ||7|| = 1: 

= (1)" = (7 . 7)" = 2 7 " . 7 + 2 7 ' . 7' (6.21) 

Now we can state the set of differential equations, that must be solved to find an analytical 
solution to the desired optimal interpolation curve: 

Proposition 31. 

Given the control points q\ , . . . , E H\ the interpolation curve 7 E C 4 (I, H\), where 7(i$) = q% 
for ti E I, will minimize ||k(/i)|| 2 g?/i if the following requirements are met: 

1. K(tl) = 

2. = 

3. For each interval ti < h < ti+i : 

7"" + (V • 7')"7 + 2(7' • 7')V + 2(7' • 7'h" = 

Now "all" that remains is to solve the above fourth-order differential equation with the given 
boundary values. This, unfortunately, is not possible with the existing mathematical knowledge. 
This is confirmed by J0rgen Sand 14 and Gerd Grubb 15 . 



6.3.6 Minimizing curvature in H\: Continuous, semi-analytical solution 

The strict mathematical derivation of the desired optimal interpolation curve stranded on an 
unsolvable differential equation. 

We can, however, again seek inspiration in the plane. The corresponding differential equation 
is here 7"" = 0, i. e. that the fourth derivative of the curve must be zero between the control 
points. This corresponds to a third-order curve (a spline). If the solution is constrained to be a 
third-order curve, the equations can be solved more easily. 

14 Associate Professor at the Institute of Computer Science, University of Copenhagen, specializing in the 
solution of systems of equations. 

16 Professor at the Institute of Mathematics, University of Copenhagen, specializing in differential equations. 
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We could therefore restrict which family of functions we would like the interpolation curve to 
belong to. This could, for example, be a kind of cubic splines. With this restriction on the set 
of possible solutions, the optimization problem could be solved. 

Depending on the choice of definition of the family of curves, this strategy for finding a solution 
will give a result corresponding to the basis for the construction of Squad. 

To keep this report at a manageable size, we will not pursue this line of thought any further. 

6.3.7 Minimizing curvature in H\\ Discretized, numerical solution 

Unfortunately, we cannot give an analytical expression for the optimal interpolation curve. 
Therefore we will try to solve a discrete version of the problem using a numerical method. 

We rewrite the problem in an equivalent discrete version. We then solve the new version of the 
problem. 

Discretization 

Given the control points Qi, ■ ■ ■ ,Qk € Hi we sought an analytical solution j(t) E C k {l,H\) 
such that there existed t\, . . . ,tk E (/), that 7(i$) = Qi, and such that the following expression 
was minimized: 

rtk 

if (7) = / \\ K (j,t)\\ 2 dt (6.22) 
Jti 

In the discrete version we will therefore attempt to solve the following problem. Given control 
points Qi, ■ ■ ■ ,Qk € H\ we seek qi, . . . , q^ E Hi such that qi t = Q t for t = 1, . . . , k and such 
that the following expression is minimized: 

N 

E = Y J KQi) \\~<Qi)\\ 2 (6-23) 
i=i 

The integral has been replaced by a sum. The "parameter width" of the interval we integrate 
over is termed It can be expressed as a centered average of the parameter width in the 

intervals immediately before and after the «'th quaternion: 

,/ \ \\<li - Qi-l\\ + \\<li - Qi+l\\ taoA\ 

H?i) = 2 ^ ' 

Another measure of the parameter distance between % and qi-\ in the approximation for l{qi) 
could be $i, where cos#j = % • %+i, instead of \\qi — qi-i\\. 

In equation 6.23 k is the discrete version of the local curvature: 




(6.25) 



64 



Note that the denominator % • % is not omitted as in the definition of local curvature (page 59). 
This is because the interpolated quaternions cannot be expected to be unit quaternions (this is 
explained below). Therefore, the denominator is not in general equal to 1; thus it cannot be 
omitted. 

In equation 6.25 the second derivative of the discrete approximation of the interpolation curve 
is used. A good centered approximation of the second derivative is ([Kincaid & Cheney, 1991], 
[Barr et al., 1992]): 

* = — ^ — (6 ' 26) 



Gradient descent 



The above equation (equation 6.23) can be minimized using gradient descent. In general terms, 
gradient descent can be described as follows. Consider the function that is to be minimized 
(commonly termed the energy function) as a hilly landscape, where the function value is the 
height of a hill at each set of coordinates. The gradient of the function in that point points 
in the steepest possible up-hill direction. Gradient descent is based on an initial guess at the 
solution. From the initial guess the gradient is computed, and a small step is taken in the 
opposite direction of the gradient (i. e. down-hill), resulting in a new point. This process is 
repeated with the new point until the function value does not become smaller by taking a step. 
In this fashion, the gradient descent method yields an approximate local minimum. It is outside 
the scope of this report to describe the theory of gradient descent any further, but we will 
describe the method in enough detail so that readers without prerequisites in the field will still 
be able to follow the derivations. 

When using numerical approximation methods such as gradient descent, it is often difficult to 
maintain restrictions on the solution space. Instead a term is added that makes the solution 
more "expensive" 16 if the solution lies outside the desired solution space. Thus, we seek a 
function that can determine if the discrete version of the interpolation curve lives in Hi, and 
thus consists of unit quaternions. This can be done by: 

g(q) =q.q-l (6.27) 

Note that Hi = {q E H | g(q) = 0}. This measure for determining if the quaternions are unit 
quaternions can be combined with the energy function E as follows: 

N 

F = Y,l(Qi) \\HQi)\\ 2 + c g(qi) 2 (6.28) 

2 = 1 

Assuming c E M suitably large, the energy function F will have a minimum approximately where 
E has a minimum, and where all % are approximately unit quaternions. 

16 We here assume "the cheaper, the better", which is not always the case in real life. Sometimes you have to 
pay extra for quality. Red wine is a good example of this, although there are, of course, exceptions in this case, 
too. 
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We want to find a minimum for F in equation 6.28 using gradient descent. We therefore have to 
find the gradient. Below, we will write q^ x for the a;'th coordinate (regarding a quaternion as a 
four-dimensional vector) in the «'th quaternion in the discrete interpolation curve. The gradient 
is 4n-dimensional. Each coordinate can be written: 



dF d 



dqi, x dqi 



N 

J2 l (<lj) \\K(<lj)\\ 2 + c g(qj) 
.i =1 



In equations 6.24, 6.25, and 6.26, % is a term in the discrete versions of 

ft((7i-i)j i^ili) an d K{Qi+i)- The corresponding terms must appear in the partial derivative: 

OF 

g— = g— \\Hm-i)\\ 2 + i{gi) WHqM 2 + Pfe+i)!! 2 + c g(qi) 2 ) 

d 

= g^— (l(Qi-i) • + %i) • + Kft+i) • + c #(%) 2 ) 

dl(<H-i)~ f \ ~ t \ i dl{qi) _ dl{q i+ i) 

2l(q i _ 1 )K(q i - 1 ) • — \-2l(q i )K(q i ) • — h 2%+i)k(%+i) • — h 

<"/;..,■ "</;..,• "</;..,• 

2c g[fu) d -i^ (6-29) 



Below we will derive the partial derivatives for the sub-expressions of the above expression. 

We introduce the notation \ x to be the vector with 1 in the a;'th coordinate, and in the other 
coordinates. We can now look forward to deriving the partial derivatives of g(qi), l{qi-i), l(Qi), 
l(Qi+i), q"-v ft"' Qi+n K(Qi-i) and k(q i+1 ): 



9 g(qj) = _d_ 
dqi, x dQi,x 

= 2qi 



(Qi • Qi ~ 1) 

d qj 
dqi, x 



— 2 q^ x 

d l(qj-i) _ d \\qj-i - gj_ g_[j + 
<9<?i,x dQi,x 2 

1 a 



(6.30) 



2 



(%-i - %-2) • (qi-i - q%-2) + \J{qi-\ - q%) • (q%-i - %)) 



lj; • (%-l - Qi) 



y/(Qi-i ~ Qi) • (Qi-i ~ Qi) t 
lx ' (qi ~ Qi-i) 



2 hi-i 



(6.31) 
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d l(qj) d \\qj - gj_ijj + \\qj - q i+ i\\ 

dqi, x dQi,x 2 

1 d 



2 <9ft, a 



(y^i -ft-i) • (ft -Qi-i) + \Jiq~i - ft+i) • (ft - ft+i)) 



1 / 1^ • (ft - gi-i) li • (qi - qi+i) 



{qi - qi-i) • {qi - ft-i) J {qi - qi+i) • (ft - ft+i) 



2 

la • (ft - ft-i) + la • (ft ~ ft+i) , g 32 . 

2 lift - ft-i|| 2 \\q i+ i -qi\\ 
d /(ft+i) 5 ||% + i - ^11 + ||% + i - q i+2 \\ 



dqi, x dq i:X 
1 d 



[yj (ft+i - ft) • (ft+i - ft) + \J (ft+i - ft+2) • (ft+i - ft+2)) 



2 <9ft v 

1 / lg • (ft+i - ft) 



2 



(ft+i - ft) • (ft+i - ft) 



lg • (ft ~ ft+i) 
2 ||ft+i - ft || 
d q"-i d - 2%_i + qi 



(6.33) 



<9ft,x <9ft,x ^(ft-i) 

'dJ~(ft-2 — ^'Ji-l t ; - v'i(-l' - -<ji-l t '/OT;^ 



Uft-i) 2 a£-(ft-2 - 2ft_i + ft) - (ft-2 - 2%_i + ft)a£-/(ft-i) 2 



Kft-i 

Kft-i) 2 a|jft - (ft-2 - 2ft-i + ft)2/(ft_i) al^^ft-i) 



Uft-i! 



4 



Uft-i) 2 !g - 2 - 2ft_i + qi)l{qi-i)-^—l{qi-i) 
l{qi-i) 4 

d q'j = d ft-i - 2gj + q i+ i 
dqi, x dqi,x l{qi) 2 

^ % ) 2 afc(*- 1 ~~ 2qi + qi +^ ~ fe- 1 ~~ 2qi + 



(6.34) 



*(ft) 4 

l{qi) 2 a§-{-2qi) ~ (ft-i - 2ft + ft+iMftJg^ft) 



/(ft) 4 

2/(%) 2 1 B + 2( ft _i - 2% + ft+i)Z(ft)^Z(ft) 
Kft? 



(6.35) 
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i+i 



dqi, x 



d k{qi) 



dqi,x 



dqi, x 



d qi- 2q i+1 + q i+2 



dqi,x 



KQi+if 



i(qi+i) 2 -^-(qi - 2<?i+i + qi+2) - {qi - 2<? i+ i + % +2 )g|^/fe+i) 2 



l(q i+ i) 4 



Kn+if ''a§-qi - (qi - 2 <?i+i + ^+2)2/(^+1) -^-Kqi+i) 



l(Qi+iY 



l{qi+i) 2 lx - 2(% - 2q i+1 + q i+2 )l{q i+ i) ^-l(q i+ i) 



l(q i+ i) 4 



a 



dqi, x 
dq'U 



-1 * 



li-l • Qi-l 



li-l 



3d 

dqi, x 



li-l 



9qi, x 
d 
dqi, x 
dq'l 
dqi, x 
dq'l 



Qi-l • Qi-l 



li-l 



Q> t 



q'l • qi 



q'i • qi 9 qi 



d (q 



qi • qi dq iyX dq iyX \ qi • q, 



9qi, x 


qi • qi 


dq'l 


Qi • H 


dqi, x 


qi • Qi 


dq'l 


Qi • Qi 


dqi, x 


Qi • qi 


dq'l 


Qi • Qi 



dqi,, 

( K'n 



dqi • qi 

dqi, x 



(ML 

\dqi,v 



(qi ' qi) 2 



\ 

! 

/ 



2 ft 



dqi 
dqi, x 



9qi,x 

d 



( Q'i' 
\dqi, x 



(6.36) 



(6.37) 



(qi ' qi) 2 

qi + q'l • ) % • qi - 2 • l^) g" • % 



(% • qi) 2 



ft(6.38) 



5 ~" 



ii+l 



li+1 



+ 1 



dqi, x 



Qi+i • %+i 



7i+i 



(6.39) 



We now have simple (but long) expressions for all the terms in the gradient. We return to 
the gradient equation (equation 6.29). All the constituent terms have been derived, and the 
gradient can be explicitly computed by substituting the derived terms for k(%_i), k(%), and 
R(qi + i) (equation 6.25) and the partial derivatives dK Q q h ~ 1 \ 9 gq 9 ^ an d dK Q q h+1 ^ (equations 6.37, 
6.38, and 6.39) into the equation of the gradient (equation 6.29). 



68 



The algorithm 



We have now derived an explicit discrete expression for the gradient, and we can perform gradient 
descent according to the following algorithm: 

Initialization 

Give a good initial guess q° = (qi, . . . , q^). An obvious way of doing this is by using one 
of the methods described earlier (Slerp or Squad). The better the initial guess, the faster 
the method converges. 

Iteration 

Write out the gradient using equation 6.29: 

VF - 9F 9F 9F 9F 9F 9F 9F 9F 

dq\,i ' <9<?i,2 ' <9<?i, 3 ' <9<?i, 4 ' '"' dq N ,i ' dq Nj2 ' dq N ^ ' dq NA 

(We set the gradient to zero in all key frames.) 

Perform the iteration on the solution guess: q 1+1 = q 1 — eVF, where e defines the step 
length 

Alternatively: Perform the iteration on the solution guess: q 1+1 = q 1 — e yy^y ■ This 
method can provide greater numerical stability. 

Termination condition 

Repeat the second step until a suitable termination condition has been reached. The 
termination condition can be dependent on the number of iterations, the total energy, F, 
the energy F in relation to the number of control points, the size of the gradient, or more 
advanced strategies. 



Alternative methods for computing the norm of the interpolation curve 

As noted above, we want the interpolation curve to lie in the space, Hi, of unit quaternions. This 
is ensured by using the penalty function g from equation 6.28. It can be difficult to determine 
when the weighting factor, c G K, of the penalty function (equation 6.28) is "suitably large." In 
theory, c must be infinite to ensure that qi £ H\. Several methods exist to handle this problem: 

Projection 

There is a simple and well-defined connection between points inside and outside the solu- 
tion space. As described previously, all quaternions on a line through the origin perform 
the same rotation. Thus the solution guess can simply be projected into the solution 
space by normalizing the generated quaternions. This can be done in each iteration or, 
alternatively, after the last iteration. 

Lagrange Multipliers 

The penalty function g{qi) can also be introduced with a Lagrange Multiplier Aj in the 
following manner: 

N 

F = E + Y, Atf(ft) 

i=l 
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A solution to equation 6.22 is a singularity for F. However, the singularity will not be a 
minimum, but instead a saddle point. Thus gradient descent is still an applicable method 
for the constituent quaternions, but gradient ascent must be performed on the auxiliary 
variables Aj. 

By including the penalty function both in equation 6.28 and a Lagrange Multiplier, the 
method become more numerically robust and converges faster. The details of this method 
can be found in [Piatt & Barr, 1988] and in [Barr et al., 1992]. 

Polar coordinates 

In the above method of solving the problem, we have reformulated the problem such that 
the restrictions on the solution space are integrated into the energy function that is to 
be minimized. However, we can restate the restriction on the solution space such that 
the complexity of the expression to be minimized does not increase. This can be done 
be representing the quaternions using polar coordinates. Thus a quaternion is written 
q = [r, (p, (9, (/>)], where r is the radius and p, 6, and cj) are the three necessary rotation 
angles. 

Using this representation, it is very easy to maintain the restriction on the solution space. 
This can be done by not performing gradient descent on the radius coordinate, that is kept 
at a constant value of 1. Thus the restriction on the solution space is maintained. 

Amongst the above representations, polar coordinates will ensure that the interpolation curve 
stays in the space, Hi, of unit quaternions. The use of Lagrange Multipliers would ensure a 
more robust and faster converging algorithm. We will not pursue either possibility, since we want 
the algorithm to be as simple as possible. We regard normalizing the generated quaternions as 
"cheating" seen from a theoretical viewpoint. We will therefore not discuss any of the three 
alternative methods any further. 

Extensions to the algorithm 

Generally, gradient descent is a method that is fairly easily expanded. The desired property 
of the interpolation curve must simply be described as a zero-crossing for some function. For 
example we might want to ensure constant angular velocity across the entire interpolation curve. 
This can be obtained using yet another penalty function: 

w(<li) = hi-i -Qi\\~ \Wi ~ Qi+i\\ (6- 4 °) 

Thus the penalty function increases with the difference between the previous and the following 
step length. This penalty can simply be introduced into the algorithm by introducing it in the 
original energy function (equation 6.28): 

N 

F = ^KQi) ||k(%)H 2 + c g 9{Qi) 2 + c w w(qi) 2 
i=i 

An extra term must be added to the gradient. This is easily derived since w{qi) 
the derivative can thus be obtained using equation 6.32: 

d w(qj) _ l x ♦ {li ~ gi-i) ]x ♦ (gj ~ Qi+i) 
dqi, x \\Qi-Qi-i\\ \\Qi+i-Qi\\ 



(6.41) 
= 21 (qi) and 

(6.42) 
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Adding the above penalty does not ensure constant speed in the interpolation curve. This is due 
to the fact that the gradient contains other terms that affect the distance between the individual 
key frames. As stated earlier, local curvature is denned such that both the geometric curvature 
and the size of the angular acceleration of the interpolation curve are included. In practice, this 
means that the angular velocity of the interpolation curve will be smaller around the key frames, 
where the curvature of the interpolation curve is maximal 17 . 

When the gradient contains terms that act in opposite directions, the algorithm becomes less 
robust. Thus setting c w "suitably large" will not necessarily ensure approximately the same 
distance between the interpolated frames. Instead this could lead to the method becoming 
numerically unstable, producing unwanted results. 

The implementation 

The above description of the algorithm is purely theoretical. Since we did not want to analyze 
the convergence and stability properties of the algorithm theoretically, there is no guarantee 
that the method works in practice. We will therefore present a number of considerations to take 
into account when implementing the algorithm. 

Multi-step minimization 

In practice, the algorithm behaves badly when given many frames. Each frame is only affected 
by its neighbors, and only key frames are positioned correctly initially. Thus the adaption to 
the key frames must propagate through all the frames between the keys and a given frame. 
The more frames that lie in-between, the more iterations are necessary for the system to attain 
an energy minimum. This give a practical upper bound on the acceptable number of frames 
between each key frame. 

An efficient way of solving this problem is by minimizing in several steps. In the first step 
relatively few frames are used, and the system attains an energy minimum after few steps. The 
frames computed in the first step are used as key frames in the second step. In the second step 
more in-between frames are added, and an energy minimum is again attained after a few steps. 
This process may be repeated an arbitrary number of times, until the desired number of frames 
has been reached. 

In figure 6.15 the result of the algorithm is seen with all frames placed during the first step. This 
gives a bad approximation curve even after many iterations. The result bears great resemblance 
to Slerp, which was used to generate the initial guess passed to the optimization algorithm. 

If the multi-step algorithm is used, a much nicer interpolation curve is produced after fewer 
iterations. In figure 6.16 the intermediate result with few frames can be seen. This is used as 
the initial configuration for the second step, and the final interpolation curve can be seen in 
figure 6.17. The angular velocity graph is also nicer when the multi-step algorithm is used. If 
the one-step method is used (see figure 6.15), the velocity curve resembles an electrocardiogram, 
while the velocity curve for the multi-step algorithm (6.17) is somewhat nicer, but still somewhat 

17 This is natural seen from a physical perspective. It is also necessary to drive slower through a sharp bend in 
the road than it is on a straight section. 
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Figure 6.15: One-step iteration. The interpolation curve contains 350 frames, and is shown 
after 500 iterations. 



uneven. This is a weakness in our implementation, but we expect that a more robust algorithm 
with better convergence properties will yield nicer velocity graphs. 

All in all, 450 iterations are used for the multi-step method versus 500 in the original version. 
The result is obviously nicer when the multi-step method is used. 
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Figure 6.16: The multi-step algorithm, first step. In the first step only 22 frames and 200 
iterations are used. 

Simplifying the parameter width - l(qi) 

In equation 6.24 a numerical approximation to the "parameter width" for the parameter of the 
interpolation curve is given. In practice it turns out that this expression has no influence on 
the shape of the interpolation curve. The algorithm becomes less numerically stable with the 
expression, however. We therefore use the simpler expression: 



(6.43) 
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Figure 6.17: The multi-step algorithm, final result. In the second step, 110 frames are in- 
terpolated using 150 iterations. In the third and final step, 550 frames and 100 iterations are 
used. 



d l(qi) _ d 



dqi, x 



dqi 







The expression for the gradient (equation 6.29) is simplified correspondingly: 



dF 

dqi, x 



2«(ft-i) • ' +2K{ qi ) . — ^l + 2k(q i+l ) - 



+2c g(qi) 



dqi, x 
d gjqi) 
dqi, x 



dqi 



dqi 



(6.44) 



(6.45) 



Weighting the curvature at the key frames 

The analytical version of the minimization of curvature ensures that the curvature is minimized 
by definition. The discrete numerical approach only gives an approximation of the solution 
with minimized curvature. The validity of the approximation to the solution depends on the 
approximations to the constituent mathematical expressions. For example the second derivative 
of the interpolation curve is approximated with: 



-i - 2qj + q i+ i 

m 2 



It is this approximation that is at the core of the approximation of the local curvature of 
the interpolation curve. Offhand, it is difficult to predict if this approximation introduces 
"weaknesses" to the numerical solution when interpolating between key frames with certain 
properties. 
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In practice the numerical method is somewhat sensitive to sharp curves. For example the curve 
in figure 6.18 is too "sharp" in the middle key frame. 

The approximation to the curvature does not adequately propagate across key frames. Let us 
reexamine the gradient (the simplified version, equation 6.45): 

- 2K( ft _i) • a yqz 11 + 2k (ft) • + 2«(ft+i) • + 2c g( qi ) ' 



oqi,x oqi,x oq iyX dq iyX dq^ x 

The last term makes sure that the quaternions remain unit quaternions. This can be ignored. 
This leaves three considerations when determining the gradient, and thus the movement of the 
individual frame during the iteration. These are the changes in curvature in the quaternion 
itself, and in both its neighbors. 

It is tempting to think that it is enough to consider curvature in the quaternion itself. However, 
if the neighbors are not taken into consideration, the result will be trivial, namely Slerp. Here 
every frame except the key frames has zero curvature. Therefore no compensation will be made 
for the large curvature in the key frames. 

Thus minimizing curvature at the key frames depends only on that the immediate neighbors of 
each key frame must "take note of" the curvature of their neighbors. This is ensured by the 
terms 2«(ft_i) • 9 K q^^ and 2k{qi + \) • 9 As shown by the example in figure 6.18, this 

is not quite sufficient. 

The approximation of the second derivative of the interpolation curve is therefore not a good 
approximation in this particular case. The approximation is too local, and should have a wider 
domain. Since this is not a report on numerical calculation methods, we have not attempted to 
find an optimal approximation expression. 

The problem we have described can be eliminated simply, though. Each key frame can simply 
"ask" its neighbors to be more "considerate." This means that we add a weighting function to 
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each expression in the gradient. The weight is dependent on whether or not a given frame is a 
neighbor of a key frame. For example, 2k(^_i) • 9 K q^'~^ will be replaced by Ck2R(c[i_i) • 9 g^'" 1 ^ , 
if qi-\ is a key frame. In practice, a good value for is about 1.2. The effect of the weighting 
function can be seen in figure 6.19. 




Figure 6.19: The multi-step algorithm with a set of key frames with a sharp curve. 200 frames 
are interpolated, and the curvature around the key frames are weighted with a factor 1.2 

Adding a constant to the algorithm "by eye" is at odds with our stated purpose of deriving 
an interpolation curve from objective criteria. The alternative is to analyze the properties of 
different numerical approximations. This, like the properties of convergence and stability, is 
outside the scope of this project. 

Remaining details 

We have not described the termination requirements for the algorithm described above. In the 
implementation we have chosen the simplest possible: The number of iterations. 

Correspondingly, we have not denned the initial guess that the algorithm uses. We have elected 
to use Slerp between each pair of key frames. This starting point is clearly not optimal. Quicker 
convergence might be achieved using Squad. Since the purpose was to derive a interpolation 
curve that is superior to Squad, it seemed to be more fair to use Slerp as the basis for the 
iteration. 

The remaining constants in the algorithm (the step size, e, in the iteration, and the penalty fac- 
tor, c) can all be calculated from robustness and convergence criteria. As previously mentioned, 
we will not describe these properties any further, and regard the constants as implementation 
details (see the introduction to the program in appendix C). 

The interpolation curve summarized 

The modest purpose of this chapter was to develop the optimal interpolation from objective, 
general criteria to the interpolation curve. 



75 



We first studied the more or less heuristic interpolation curves that are to be found in the 
literature. These included the naive LinMat, LinEuler, Lerp, the simple Slerp, and finally the 
convincing Squad. 

Using these simple interpolation curves as a basis we tried to define which class of functions the 
interpolation curve should belong to (page 56). We were not able to determine which definition 
of smoothness was suitable to define the class of desired functions. 

We then attempted to define an interpolation curve that minimized the integral of the local 
curvature (denned in equation 6.17) of the interpolation curve. These derivations required 
that the interpolation curve is four times differentiable with continuous derivatives (i. e. 7(f) € 
C 4 (I, H\)). This is noted on page 62 in section 6.3.5. Unfortunately, this derivation gave rise 
to a fourth-order differential equation that we were unable to solve. Thus it is irrelevant to 
consider the open questions from section 6.3.2: How many times differentiable should the curve 
be and are singularities allowed? 

Thus we settled for a discrete, numerical solution. We have presented a method based on 
gradient descent. We examined and refined the method. The final result were some very pleasing 
interpolation curves. 

As our final interpolation algorithm we will choose the basic algorithm from section 6.3.7 with 
the relative distribution of frames in the sub-intervals (see section 6.2.1 on page 55) and with 
the weighting of the curvature around the key frames. This interpolation curve we will name 
Spring (for Spherical Interpolation using Numerical Gradient descent). In figure 6.20, the 
effect of using Spring can be seen. The result is only marginally better than the version that 
does not include special treatment of curvature at the key frames. In the next chapter we will 
see examples where Spring much more clearly demonstrates its value. 




Figure 6.20: The effect of Spring. 550 frames have been interpolated and curvature is weighted 
with a factor 1.3 
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Chapter 7 

Squad and Spring 



In chapter 6 we treated a series of interpolation algorithms. The most convincing among the 
known algorithms was Squad. In this chapter we will compare Squad with our own algorithm 
Spring. 

The comparison will be based on a number of illustrative examples. 



7.1 Example: A semi circle 



First we check if Squad and Spring can produce nice rounded curves. We have placed the key 
frames as corners in a spherical square. Around the center key frames the interpolation curve 
should approximately be a semi circle. Figure 7.1 and 7.2 show that both curves meet this 
requirement nicely. 




Figure 7.1: A simple curve with soft rounded corners interpolated with Squad. A total of 550 
frames have been used in the interpolation. 
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Figure 7.2: A simple curve with soft rounded corners interpolated with Spring. A total of 430 
frames and 700 iterations in the steps have been used. No extra weight on the curvature at the 
key frames has been added (see section 6.3.7 page 73). 



7.2 Example: A nice soft curve 

This next example should be no real challenge for either of the algorithms. The expected 
interpolation curve is simply a nice rounded curve with no sharp corners. The intersection in 
the curve should pose no problem. 
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Figure 7.3: Interpolation with Squad of a soft curve with an intersection. A total of 550 frames 
have been used in the interpolation. 



Both interpolation curves on figure 7.3 and figure 7.4 are nice. However, the curve has more 
rounded corners for Spring than for Squad (even though no extra weight has been added to the 
curvature at the key frames - see section 6.3.7 page 73). This means that the curve for Spring 
has smaller curvature. Furthermore the velocity graph for Spring is more constant than for 
Squad. This implies that Spring behaves as desired. 
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Figure 7.4: Interpolation with Spring of a soft curve with an intersection. A total of 550 
frames and 700 iterations distributed over the interpolation steps. No extra weight is added to 
the curvature around the key frames (see section 6.3.7 page 73) 



7.3 Example: Interpolation curve with cusp 



Squad can produce interpolation curves with quite pointy corners at the key frames. 




Figure 7.5: A curve with a cusp interpolated with Squad. 200 frames have been used. 



The interpolation curve for Squad (figure 7.5) reveals a nasty pointy curve. However, Spring is 
able to produce a nice smooth rounded interpolation curve (figure 7.6). 

That the interpolation curve for Squad has a sharp corner does not contradict the proven fact 
that Squad is differentiable (section 6.2.1 page 51). At the key frame with the sharp corner the 
velocity of the curve is zero. Thereby the function Squad remains differentiable although the 
geometric appearance of the curve is not intuitively differentiable. 
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Figure 7.6: A curve with a cusp interpolated with Spring. A total of 200 frames and 700 
iterations distributed over three steps. The relative weight of the curvature around the key frames 
is 1.3 (see section 6.3.7 page 73). 



7.4 Example: A pendulum 



We continue to investigate curves with large curvature. In this example we use a curve with 
infinite curvature - a pendulum motion. This is achieved with three key frames where the first 
and last are equal. 




Figure 7.7: Squad displaying pendulum motion over 250 frames. 



The desired behaviour of the interpolation curve is not intuitively obvious. Since all key frames 
are on a arc one would expect the curve to remain on this arc. 

The figures (7.7 and 7.8) shows the same behaviour for both curves - the pendulum motion. 
Note that Lerp and Slerp would produce the same curve. 
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Figure 7.8: Spring displaying pendulum motion over 200 frames and 700 iterations in three 
steps. The relative weight of the curvature around the key frames is 1.5 (see section 6.3.7 page 
73). 

7.5 Example: A perturbed pendulum 

Even though it is very reasonable that the interpolation curve remains on the same arc when 
the key frames are all on an arc, this is not necessarily correct. In principle the curve has 
infinite curvature at the center key frame. Since Spring is supposed to minimize curvature this 
is somewhat disappointing. 

To understand this, it is necessary to bear the algorithm in mind. Using gradient descent, each 
frame will move slightly in each step to decrease the curvature. But in which direction should 
the frames close to the center key frame move? Since the curve is symmetric, the gradient at 
each frame will be zero. The problem is that this does not imply a local minimum but a local 
maximum instead. 

Thus, the pendulum is an example of how it is possible to confuse Spring. We investigate this 
further by perturbing the first and last key frames slightly so the curve is no longer a pure 
pendulum. It is only just possible to see in the curve for Squad (figure 7.9). 

Figure 7.10 shows how the pure arc is no longer a (repelling) fix point for the gradient descent 
algorithm Spring. The minimal perturbation allows the interpolation curve to be nice and 
rounded at the center key frame. 

7.6 Example: Global properties 

The final example demonstrates a fundamental difference between Squad and Spring. In each 
interval the interpolation curve for Squad is denned exclusively from the two previous and the 
two following key frames — i. e. a local definition. In contrast, the interpolation curve for Spring 
is globally defined. 
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Figure 7.9: The perturbed pendulum interpolated by Squad using 250 frames. 




Figure 7.10: The perturbed pendulum interpolated by Spring using 200 frames and 700 iterations 
in three steps. The relative weight of the curvature around the key frames is 1.5 (see section 
6.3.7 page 73) 

Figure 7.11 shows the interpolation curve for Squad on a set of five key frames. The first three 
key frames lie approximately on an arc and therefore the interpolation curve is an arc in the 
first interval. Likewise the interpolation curve form an arc in the last interval. 

In contrast the interpolation curve for Spring (figure 7.12) is nice and smooth. The global 
structure of the algorithm allows the curve to distribute the curvature evenly across all the 
intervals. Instead of having excessive curvature at the center key frame, a part of the curvature 
is propagated to the outer intervals. 

It should be noted that it is necessary to add a relatively large weight to the curvature around 
the key frames in this example. However, we have no doubt that an algorithm with a better 
numerical approximation for the constituent expressions (in particular q" in equation 6.26) would 
produce the same result — without having to fit the parameters of the program to the example. 



82 




Figure 7.11: Squad producing pointy curve using 550 frames. 




Figure 7.12: Spring avoiding the pointy curve using 430 frames and 900 iterations. The relative 
weight of the curvature around the key frames is 4 (see section 6.3.7 page 73). 
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7.7 Conclusion 

The fundamental differences between the to methods are displayed below. 



Property for the 


Characteristics for Squad 
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must be fitted to some examples. 



The choice between Squad and Spring is not obvious. If a simple algorithm yielding nice results in 
most cases is needed, then the simple Squad will suffice. If really nice interpolation is mandatory 
in all cases the more complex Spring will be more appropriate. 
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Chapter 8 

The Big Picture 



In this final chapter we will first attempt to discuss our work in relation to the available literature 
on quaternions and interpolation of rotations. Using this as a starting point we will point out 
relevant topics for future work. 



8.1 Comparison to previous work 

This paper has covered five main topics: 

• Rotation modalities (Section 3). 

• Quaternion mathematics (Section 3.3). 

• Curves for interpolation of rotations — Heuristic approach (Section 6.2). 

• Curves for interpolation of rotations — Analytic minimization of the local curvature (Sec- 
tion 6.3.5). 

• Curves for interpolation of rotations - Numerical minimization of the tangential curvature 
(Section 6.3.7). 

For each of the main topics we will summarize our contributions compared to the previous work. 
Quaternion mathematics 

Quaternion mathematics has been treated several times in the literature ([Hamilton, 1853], 
[Hamilton, 1899], [Pervin & Webb, 1992], [Shoemake, 1994b], [Maillot, 1990], and [Kim et al., 1996]). 
[Pervin & Webb, 1992] should be noted for a treatment of the basic mathematical properties — 
including the logarithm and exponential functions. In [Kim et al., 1996], a general framework 
for differentiation is given (though the reader is referred to a differential geometry text for the 
proof), and the derivative of exp is derived based on this framework. None of the articles have 
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given a complete treatment of the necessary mathematics. Differentiation of the quaternion 
functions is very central in the study of the smoothness of the interpolation curves. 

This report includes a comprehensive treatment of quaternion math. In particular, we have 
derived all the differentiation equations necessary for proving the desired properties of the in- 
terpolation curves (in section 3.3.8). 

Curves for interpolation of rotations — Heuristic approach 

Interpolation curves in the plane have inspired several interpolation curves for rotations. The two 
most important are Slerp [Shoemake, 1985] and Squad [Shoemake, 1987]. Shoemake attempts 
to prove the differentiability of Squad but the proof is flawed 1 . In [Kim et al., 1996], a more 
general result is given that entails the differentiability of Squad. Our result is derived using less 
advanced mathematics, and may be more easily accessible. 

This is the first report which contains a comprehensive treatment of the two most important 
heuristic quaternion curves (section 6.2). All the known expressions for Slerp are stated and 
the correctness of their properties is proven. For Squad we apply the derived differentiation 
equations to prove the differentiability of the function. 

Other heuristic approaches 

Like Squad, a number of quaternion interpolation curves have been made from general spher- 
ical cubic curves. These curves are fitted to the control points, thus yielding relatively nice 
interpolation curves. Most often these curves are inspired by cubic curves in the plane. 

Examples of this are the fairly simple spherical Catmull-Rom B-spline in [Schlag, 1994] and a 
spherical Bezier curve in [Shoemake, 1985]. 

In this report, we have not investigated this approach. 

Curves for interpolation of rotations — Analytic approach 

The first attempt to derive an optimal interpolation curve from a set of objective criteria can be 
seen in [Barr et al., 1992]. The paper states an expression minimizing the tangential curvature 
(in this paper also called the local curvature). However, the paper makes no attempt to state or 
solve the differential equation, that corresponds to the optimal curve. 

In this paper we state a set of objective criteria for the optimal interpolation curve (section 
6.3.5) — inspired by [Barr et al., 1992]. We then derive the fourth order differential equation 
whose solution will minimize the local curvature. Unfortunately, we are not able to solve the 
differential equation analytically. 



1 As mentioned earlier [Shoemake, 1987] is unavailable but Shoemake himself has stated [Shoemake, 1997] that 
the proof was flawed. 



86 



Minimization of the local curvature — Numerical approach 

A brief description of a numerical solution to the problem of minimizing the tangential curvature 
is sketched in [Barr et al., 1992]. [Piatt & Barr, 1988] contains a more sophisticated method 
(using Augmented Lagrangian constraints in a Finite elements method). A method yeilding 
faster convergence is presented in [Ramamoorthi & Barr, 1997]. 

We give a complete algorithm for finding a discrete solution to the numerical version of the 
minimization problem (section 6.3.7). 

Complete treatment 

To our knowledge, there exists no other complete treatment of quaternion mathematics and the 
applications in interpolation of rotations. 

This report combines a comprehensive overview including both a thorough treatment of the basic 
quaternion mathematics and as well, the most important methods for interpolating orientations 
is space. This is where quaternions really show their strength (see sections 3.3.6, 3.4, and chapter 
6). 

8.2 Future work 

Obviously, it would be very nice to derive an analytic solution of the differential equation that 
minimizes the local curvature of the interpolation curve. Since differential equations are centuries 
old this is not very likely to happen in the immediate future. Therefore it would be more realistic 
to establish a more robust numerical solution with a faster convergence. This is beyond the scope 
of this report. An excellent starting point for this approach would be [Piatt & Barr, 1988] and 
[Ramamoorthi & Barr, 1997]. 

Another direction could be the development of more specialized applications. As an example, 
the movements of a camera should not necessarily be interpolated in the same manner as the 
moving object during animation. A stable horizon is possibly a desired feature for the camera 
(i. e. the camera must not tilt upside down). Relevant introductions to this line of work are 
[Shoemake, 1994b] and [Shoemake, 1994a]. 

Another example of specialization is applications where the interpolation need not be smooth. 
For instance, in the plane the interpolation of a fly or a UFO need not be smooth either. An 
example of extraction of certain properties of interpolation curves in 3D (tension, continuity 
and bias control) can be found in [Kochanek & Bartels, 1984]. 

Finally, yet another example of specialization of the interpolation curve can be studied in 
[Barr et al., 1992]. In this method, the angular velocity can be given explicitly at the first and 
last key frame (the authors call this Angular Velocity Constraints). An obvious generalization 
of this would be the ability to supply the angular velocity at an arbitrary key frame. 
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Appendix A 



Conventions 



In this report we have used following conventions: 



Coordinate system 

We use a right-handed coordinate system. In computer graphics it is common to use 
a left-handed coordinate system. This allows the z-axis to point "into" the screen which 
seems natural. Since we primarily use coordinates for mathematical derivations we have 
chosen to use the mathematical standard — the right-handed coordinate system. 

Rotation 

Still using the mathematical standard we rotate counter-clockwise. The direction of ro- 
tation about an axis is obtained by the right-hand rule: Hold the axis with right hand 
and the thumb pointing in the positive direction of the axis. A positive rotation will now 
rotate in the direction of the fingers (apart from the thumb). 

This is illustrated below: 



y 




Rotation about z brings x into y 
Rotation about y brings z into x 
Rotation about x brings y into z 



z. 



Euler angles 

Rotation by Euler angles is defined by a rotation about each of the three coordinate axes. 
To make this unambiguous it is necessary to define the order of rotation. 

The specific order of rotation is of no importance in this paper and we arbitrarily choose 
x, y, z. Other conventions are described in [Craig, 1986]. 
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Appendix B 



Conversions 



In this appendix we show the conversions between different representations for rotation: Euler 
angles, matrices, and quaternions. 



B.l Euler angles to matrix 



Rotation about the £-axis by the angle a followed by rotation about the y-axis by the angle f3 
concluded by rotation about the £-axis by the angle 7 is written in matrix 1 notation: 



cos 7 
sin 7 







— S1117 
cos 7 





cos/3 


— sin (3 




sin/3 


cos/3 






cos a 
sin a 






— sin a 
cos a 




cos /3 cos 7 cos 7 sin a sin /3 — cos a sin 7 

cos /3 sin 7 cos a cos 7 + sin a sin (3 sin 7 
— sin /3 cos /3 sin a 





cos a cos 7 sin /3 + sin a sin 7 
- (cos 7 sin a) + cos a sin /3 sin 7 
cos a cos /3 




B.2 Matrix to Euler angles 



The rotation matrix derived above is the starting point for the conversion from matrix to Euler 
angles. The conversion to Euler angles requires the inverse trigonometric functions. Neither 
arcsin nor arccos will, by itself, yield values on the entire interval from — tt to tt. We therefore 
want to establish both the sin and cos values for all the angles. Using both the sin and cos 
values the original angles can be determined in the interval ] — tt, tt] as v = sgn(sinu) arccos (v), 
where sgn is defined: sgn(O) = and sgn(a;) = jjj^. 

1 Using homogeneous matrices the rotation matrices are 4x4. 
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We can directly determine sin/3 as —Ru- Isolating cos/3 yields the following equations: 



COS 2 /3 = (-R U R 32 R 31 -R 33 R 21 )/R 12 

COS 2 /3 = (-i? n i?33#31 +i? 32 i?2l)/i?13 

COS 2 /3 = {-R 32 R n R 21 +R U R33)/R22 

COS 2 /3 = {-R 31 R^R 21 -R U R32)/R23 

From this it is not possible to determine the sign of cos /3. This means that we might as well 
determine /3 directly from sinf3. Either way it is only possible to determine /3 in the interval 
[— f , f] corresponding to the assumption that cos/3 is positive. 

Determining cos and sin for a and 7 requires cos /3. If the assumption that cos /3 is positive does 
not hold, then also a and 7 are determined incorrectly. Unfortunately this is the best that can 
be done. 

The equations for cos and sin for each of the constituent angles are shown below. From this the 
angles can be determined as stated above. 



cos a 
sin a 
cos 7 
sin 7 



/3 = arcsin(— i?3i) (Assuming /3 € 
R33 



IT IT 

2' 2 



cos /3 
R32 
cos /3 
Rn 
cos /3 
R21 
cos /3 



Obviously this requires that cos/3 7^ 0. If cos/3 = then /3 = ±|. This corresponds to gimbal 
lock (see section 4.1) and therefore it is impossible to distinguish a from 7. Thus we arbitrarily 
define 7 = 0. In this situation the angles can be determined: 



/3 = arcsin(— i?3i) (assuming/3 € 

cos a = R 22 

sin a = —R 2 3 

7 = 



7T 7T 

2' 2 



B.3 Quaternion to matrix 



Rotation of the vector p = (x, y, z) with the quaternion q is done by the operation q [0,p] q -1 . 
We want to determine the corresponding matrix which multiplied on [x y z 1] T from the left 
will yield the same result. 
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The product of two quaternions q v = [w, (a, b, c)] and = [s, (x, y, z)] (written in i, j, k-notation) 
is: 

Qvqii = {w + ia + j& + kc)(s + \x + jy + kz) 

= (ws — ax — by — cz) + i(as + wx — cy + bz) + 
j(bs + ivy + cx — az) + k(cs + wz — bx + ay) 

Written as columns using sloppy notation 2 this equals: 



QvQh = 



a 




X 




wx — cy + bz + as 


b 




y 




cx + wy — az + bs 


c 




z 




—bx + ay + wz + cs 


w 




s 




—ax — by — cz + ws 



From this we can write the matrices corresponding to multiplying from the left and from the 
right with a quaternion. First we determine V qv such that V Qv qh = q v Qhi where q v and q v qh are 
written as columns: 



V, 



qv 



w —c b a 

c w —a b 

—b a w c 

-a —b —c w 



Then we write Hq h , such that Hq,q v — qvQh- 



H, 



Qh 



s 


z 


-y 


X 


—z 


s 


X 


y 


y 


—x 


s 


z 


—x 


-y 


—z 


s 



We are now ready to write the matrix M, such that Mp 
and q~ l = [s, (—x, —y, —z)] we get: 



q [0,p] q 1 . Using q = [s,(x,y,z)] 



M 



V q H q -i 



s 


— z 


y 


X 




s - 


-z 


y -x 


z 


s 


—x 


y 




z 


s 


-x -y 


-y 


X 


s 


z 




-y 


X 


s —z 


—x 


-y 


—z 


s 




X 


y 


z s 


1 - 


2(2/ 2 


+ z 2 ) 


2xy — 2sz 




2sy + 



2xy + 2sz 
—2sy + 2xz 




-2(x 2 + z 2 ) 
2sx + 2yz 




-2sx + 2yz 
-2(x 2 + y 2 ) 







n 
l 



2 The quaternion [s, (x,y,z)\ is written as the column [x, y, z, s] T , and x 9 denotes quaternion multiplication 
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B.4 Matrix to Quaternion 



Conversion from a rotation matrix to the corresponding unit quaternion uses the matrix M 
derived above. First we find s: 



Mn + M 22 + M 33 + M 44 = 4 - 4(x 2 + y 2 + z 2 ' 



4-4(l-s 2 ) Da s 2 +x 2 + y 2 + z 2 = 1 
As 2 



This yields s . Now the other values follow: 

S = ± i + M 22 + M33 + M 

M 32 - M 23 



a: 



4s 

M13 - M31 
y = 4, 

M 2 i - M12 

Z = 4, 



The sign of s cannot be determined. Depending on the choice of sign for s the signs for x, 
y and z change as well. This means choosing between a quaternion and the corresponding 
negative quaternion. These quaternions yield the same rotation but the interpolation curve can 
be influenced by this choice. 

Since the entire interpolation is calculated on quaternions this will pose no practical problem. 
Therefore we simply choose the positive square root. 



B.5 Between quaternions and Euler angles 

These conversions can simply be achieved by going via matrices using the conversions stated 
above. 

Therefore the limitation on the f3 angle from the conversion between matrix and Euler angles 
will hold for the conversion from quaternions to Euler angles as well. It is not possible to avoid 
this limitation (or a similar one) using a direct conversion between quaternions and Euler angles. 
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Appendix C 

Implementation 



This chapter contains a brief description of the program we have developed for visualization 
including the standard packages we have used. 

quat displays in a window on the screen how an object (the letter "R" ) looks when it is rotated 
in space. The desired key orientations are supplied through a resource file together with various 
parameters. The interpolation method is supplied as a command line parameter. The available 
methods are: 

lineuler Linear interpolation between Euler angles (section 6.1.1). 

linmat Linear interpolation between rotation matrices (section 6.1.2). 

lerp Linear interpolation between quaternions (section 6.1.3) 

slerp Spherical linear interpolation between quaternions (section 6.1.5). 

squad Spherical spline interpolation between quaternions (section 6.2.1) 

slerpsvupti Minimization of the tangential curvature using a gradient descent 

method with Slerp as initial solution (section 6.3.7). 
squadsvupti Minimization of the tangential curvature using a gradient descent 

method with Squad as initial solution (section 6.3.7). 
justdoit Minimization of the tangential curvature using a gradient descent 

method applied three times with Slerp as initial solution (section 6.3.7). 

The curious names slerpsvupti, squadsvupti and justdoit are used for historical reasons. 

We use the graphics library SPHIGS [Sklar & Brown, 1993] for displaying the animation on the 
screen. However, we have added the ability to export the animation as a series of PPM files. The 
PPM files are converted to an animated GIF using the program convert [ImageMagic, 1997]. 
A few examples of this can be seen at http://kantine.diku.dk/~myth/gif 

The program quat produces a visualization of the interpolated quaternions and an approxi- 
mation of the velocity (see chapter 5). The velocity is displayed as a two dimensional graph 
generated by the program gnuplot [Williams & Kelley, 1993]. The visualization of the interpo- 
lation curve is made using the ray tracer POV-ray [POV, 1997]. 
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C.l The basic structure of quat 



We have used C++ for writing quat. The object oriented language allows us to implement 
classes for each mathematical concept: matrix, vector, quaternion and so on. Apart from this 
separate objects handle interpolation and visualization. 

The source code can be obtained from the authors. 
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